ANTES DE ESTUDIAR PCA

Los datos que se usan en el analisis de componentes principales en R y Python se pueden descargar en

https://github.com/jaimeisaac2020/Python-analsisis-basicos/blob/mi-github/nadadores.xlsx

Tambien puedes estudiar el PCA en R de forma muy detallada en el siguiente Rpubs: https://rpubs.com/Joaquin_AR/287787

Algo para introducir el PCA

# Definir la matriz A de 4x4
A <- matrix(c(1, 2, -1, 1, 0, 1, 3, 1, 0, 0, 2, 0, 0, 0, 1, -1), 
            nrow = 4, ncol = 4, byrow = TRUE)

# Mostrar la matriz A
A

     [,1] [,2] [,3] [,4]
[1,]    1    2   -1    1
[2,]    0    1    3    1
[3,]    0    0    2    0
[4,]    0    0    1   -1

# Calcular los autovalores de A
print(eigen(A)$values)

[1]  2  1  1 -1

## Comparar los siguientes cálculos:
sum(diag(A)) # Calcula la traza de A

[1] 3

sum(eigen(A)$values) # Calcula la suma de los autovalores de A

[1] 3

## Comparar los siguientes cálculos:
det(A) # Calcula el determinante de A

[1] -2

prod(eigen(A)$values) # Calcula el producto de los autovalores de A

[1] -2

# Calcular la transpuesta de A
t(A)

     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    2    1    0    0
[3,]   -1    3    2    1
[4,]    1    1    0   -1

# Observar que las trazas de una matriz y su transpuesta son iguales
sum(diag(t(A)))

[1] 3

# Observar que los determinantes de una matriz y su transpuesta son iguales
det(t(A))

[1] -2

# Observar que los autovalores de una matriz y su transpuesta son los mismos
eigen(t(A))$values

[1]  2 -1  1  1

# Calcular la inversa de A
solve(A)

     [,1] [,2] [,3] [,4]
[1,]    1   -2  4.0   -1
[2,]    0    1 -2.0    1
[3,]    0    0  0.5    0
[4,]    0    0  0.5   -1

# Verificar que A y su inversa son efectivamente inversas (debe dar la matriz identidad)
A %*% solve(A)

     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    1    0    0
[3,]    0    0    1    0
[4,]    0    0    0    1

# Observar que los determinantes de una matriz y su inversa son inversos
det(solve(A))

[1] -0.5

# Observar que los autovalores de una matriz y su inversa son inversos
1 / eigen(A)$values # Autovalores de la inversa

[1]  0.5  1.0  1.0 -1.0

eigen(solve(A))$values

[1]  1.0  1.0 -1.0  0.5

eigen (A) $ vectors # Calcula los autovectores de A

           [,1] [,2]          [,3]       [,4]
[1,] 0.86402765    1 -1.000000e+00  0.0000000
[2,] 0.48001536    0  1.110223e-16 -0.4472136
[3,] 0.14400461    0  0.000000e+00  0.0000000
[4,] 0.04800154    0  0.000000e+00  0.8944272

eigen (A) $ vectors [ , 1 ] # Muestra e l primer autovector

[1] 0.86402765 0.48001536 0.14400461 0.04800154

## Ver i f iquemos que es autovector de autovalor 2:
A%*%eigen (A) $ vectors [ , 1 ]

           [,1]
[1,] 1.72805530
[2,] 0.96003072
[3,] 0.28800922
[4,] 0.09600307

2*eigen (A) $ vectors [ , 1 ]

[1] 1.72805530 0.96003072 0.28800922 0.09600307

sqrt(sum( eigen (A) $ vectors [ , 1 ] ^ 2 ) ) # Calcula la norma del primer autovector

[1] 1

import numpy as np

# Definir la matriz A de 4x4
A = np.array([[1, 2, -1, 1],
              [0, 1, 3, 1],
              [0, 0, 2, 0],
              [0, 0, 1, -1]])

# Mostrar la matriz A
print("Matriz A:")

Matriz A:

print(A)

[[ 1  2 -1  1]
 [ 0  1  3  1]
 [ 0  0  2  0]
 [ 0  0  1 -1]]

# Calcular los autovalores de A
eigenvalues_A = np.linalg.eigvals(A)
print("\nAutovalores de A:")


Autovalores de A:

print(eigenvalues_A)

[ 1.  1. -1.  2.]

## Comparar los siguientes cálculos:
# Traza de A (suma de los elementos de la diagonal principal)
trace_A = np.trace(A)
print("\nTraza de A:")


Traza de A:

print(trace_A)

print(sum(np.diag(A)))

# Suma de los autovalores de A
#sum_eigenvalues_A = np.sum(eigenvalues_A)
#print("Suma de los autovalores de A:")
#print(sum_eigenvalues_A)

## Comparar los siguientes cálculos:
# Determinante de A
det_A = np.linalg.det(A)
print("\nDeterminante de A:")


Determinante de A:

print(det_A)

-2.0

# Producto de los autovalores de A
prod_eigenvalues_A = np.prod(eigenvalues_A)
print("Producto de los autovalores de A:")

Producto de los autovalores de A:

print(prod_eigenvalues_A)

-2.0

# Transpuesta de A
A_transpose = A.T
print("\nTranspuesta de A:")


Transpuesta de A:

print(A_transpose)

[[ 1  0  0  0]
 [ 2  1  0  0]
 [-1  3  2  1]
 [ 1  1  0 -1]]

# Observar que las trazas de una matriz y su transpuesta son iguales
trace_A_transpose = np.trace(A_transpose)
print("\nTraza de la transpuesta de A:")


Traza de la transpuesta de A:

print(trace_A_transpose)

print(np.trace(A))

# Observar que los determinantes de una matriz y su transpuesta son iguales
det_A_transpose = np.linalg.det(A_transpose)
print("Determinante de la transpuesta de A:")

Determinante de la transpuesta de A:

print(det_A_transpose)

-2.0

print("Determinante  de A:")

Determinante  de A:

print(np.linalg.det(A))

-2.0

# Observar que los autovalores de una matriz y su transpuesta son los mismos
eigenvalues_A_transpose = np.linalg.eigvals(A_transpose)
print("\nAutovalores de la transpuesta de A:")


Autovalores de la transpuesta de A:

print(eigenvalues_A_transpose)

[ 2. -1.  1.  1.]

print("\nAutovalores  de A:")


Autovalores  de A:

print(np.linalg.eigvals(A))

[ 1.  1. -1.  2.]

# Calcular la inversa de A
A_inverse = np.linalg.inv(A)
print("\nInversa de A:")


Inversa de A:

print(A_inverse)

[[ 1.  -2.   4.  -1. ]
 [ 0.   1.  -2.   1. ]
 [ 0.   0.   0.5  0. ]
 [-0.  -0.   0.5 -1. ]]

# Verificar que A y su inversa son efectivamente inversas (debe dar la matriz identidad)
identity_check = np.dot(A, A_inverse)
print("\nVerificación de que A y su inversa son inversas (matriz identidad):")


Verificación de que A y su inversa son inversas (matriz identidad):

print(identity_check)

[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]]

A_inverse @ A

array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

# Observar que los determinantes de una matriz y su inversa son inversos
det_A_inverse = np.linalg.det(A_inverse)
print("\nDeterminante de la inversa de A:")


Determinante de la inversa de A:

print(det_A_inverse)

-0.5

print("\nDeterminante  de A:")


Determinante  de A:

print(np.linalg.det(A))

-2.0

# Observar que los autovalores de una matriz y su inversa son inversos
eigenvalues_A_inverse = np.linalg.eigvals(A_inverse)
print("\nAutovalores de la inversa de A:")


Autovalores de la inversa de A:

print(eigenvalues_A_inverse)

[ 1.   1.  -1.   0.5]

print("Recíprocos de los autovalores de A:")

Recíprocos de los autovalores de A:

print((eigenvalues_A)**(-1))

[ 1.   1.  -1.   0.5]

print("\nAutovalores  de A:")


Autovalores  de A:

print(np.linalg.eigvals(A))

[ 1.  1. -1.  2.]

# Calcular los autovalores y autovectores de A
eigenvalues_A, eigenvectors_A = np.linalg.eig(A)

# Mostrar los autovectores de A
print("Autovectores de A (columnas):")

Autovectores de A (columnas):

print(eigenvectors_A)

[[ 1.00000000e+00 -1.00000000e+00  0.00000000e+00  8.64027649e-01]
 [ 0.00000000e+00  1.11022302e-16 -4.47213595e-01  4.80015361e-01]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  1.44004608e-01]
 [ 0.00000000e+00  0.00000000e+00  8.94427191e-01  4.80015361e-02]]

# Extraer el primer autovector (primera columna de eigenvectors_A)
first_eigenvector = eigenvectors_A[:, 0]
print("\nPrimer autovector:")


Primer autovector:

print(first_eigenvector)

[1. 0. 0. 0.]

## Verificar que el primer autovector es un autovector asociado al primer autovalor
# Multiplicar A por el primer autovector
A_times_first_eigenvector = np.dot(A, first_eigenvector)
print("\nA multiplicado por el primer autovector:")


A multiplicado por el primer autovector:

print(A_times_first_eigenvector)

[1. 0. 0. 0.]

# Multiplicar el primer autovector por su autovalor correspondiente
lambda_1 = eigenvalues_A[0]  # Primer autovalor
lambda_times_first_eigenvector = lambda_1 * first_eigenvector
print("Autovalor multiplicado por el primer autovector:")

Autovalor multiplicado por el primer autovector:

print(lambda_times_first_eigenvector)

[1. 0. 0. 0.]

## Calcular la norma del primer autovector
norm_first_eigenvector = np.sqrt(np.sum(first_eigenvector**2))
print("\nNorma del primer autovector:")


Norma del primer autovector:

print(norm_first_eigenvector)

1.0

Supongamos que deseamos explorar en nuestra población los factores de riesgo de sufrir una enfermedad coronaria. De estudios anteriores sabemos que se consideran como factores de riesgo para la enfermedad coronaria: la hipertensión arterial, la edad, la obesidad, el tiempo de antigüedad en el diagnóstico de hipertensión, el pulso, y el stress.

Para la investigación se seleccionan al azar 20 pacientes hipertensos de la población objetivo sobre los cuales se miden las siguientes variables:

*\(X_1\) : presión arterial media en \(\mathrm{mm} / \mathrm{Hg}\),

*\(X_2\) : edad en años

*\(X_3\) : peso en kg

* \(X_4\) : superficie corporal en \(\mathrm{m}^2\),

* \(X_5\) : tiempo transcurrido desde el diagnóstico de hipertensión en años,

X6: pulsaciones por minuto,
X7: medida asociada al stress.

# Cargar paquetes necesarios
library(ggplot2) # Paquete para confeccionar dibujos
library(devtools) # Colección de herramientas de desarrollo para paquetes

Warning: package 'devtools' was built under R version 4.4.3

Cargando paquete requerido: usethis

Warning: package 'usethis' was built under R version 4.4.3

# Instalar el paquete "ggbiplot" desde GitHub (si no está instalado)
#if (!requireNamespace("ggbiplot", quietly = TRUE)) {
#  install_github("vqv/ggbiplot") # Instala el paquete desde GitHub
#}

library(ggbiplot) # Paquete para visualización de componentes principales

Cargando paquete requerido: plyr

Cargando paquete requerido: scales

Cargando paquete requerido: grid

library(readxl) # Permite leer archivos xlsx

PCA EN R

# Importar la base de datos con la cual se va a trabajar
nad <- read_excel("nadadores.xlsx") # Asegúrate de especificar la ruta correcta al archivo
print(nad)

# A tibble: 14 × 5
   Nadador   tr1   tr2   tr3   tr4
     <dbl> <dbl> <dbl> <dbl> <dbl>
 1       1    10    10    13    12
 2       2    12    12    14    15
 3       3    11    10    14    13
 4       4     9     9    11    11
 5       5     8     8     9     8
 6       6     8     9    10     9
 7       7    10    10     8     9
 8       8    11    12    10     9
 9       9    14    13    11    11
10      10    12    12    12    10
11      11    13    13    11    11
12      12    14    15    14    13
13      13    10    10    12    13
14      14    15    14    13    14

# Crear un data frame con las columnas relevantes
nadadores <- data.frame(nad[, 2:5])
print(nadadores)

   tr1 tr2 tr3 tr4
1   10  10  13  12
2   12  12  14  15
3   11  10  14  13
4    9   9  11  11
5    8   8   9   8
6    8   9  10   9
7   10  10   8   9
8   11  12  10   9
9   14  13  11  11
10  12  12  12  10
11  13  13  11  11
12  14  15  14  13
13  10  10  12  13
14  15  14  13  14

# Realizar el análisis de componentes principales sin estandarización
nad.pca.cov <- prcomp(nadadores, center = TRUE, scale. = FALSE)

# Realizar el análisis de componentes principales con estandarización
nad.pca.cor <- prcomp(nadadores, center = TRUE, scale. = TRUE)

# Resumen de las varianzas explicadas por las componentes principales
summary(nad.pca.cor) # Con estandarización

Importance of components:
                          PC1    PC2     PC3     PC4
Standard deviation     1.7098 0.9573 0.34831 0.19690
Proportion of Variance 0.7309 0.2291 0.03033 0.00969
Cumulative Proportion  0.7309 0.9600 0.99031 1.00000

summary(nad.pca.cov) # Sin estandarización

Importance of components:
                          PC1    PC2    PC3     PC4
Standard deviation     3.5840 1.9858 0.7030 0.42241
Proportion of Variance 0.7356 0.2258 0.0283 0.01022
Cumulative Proportion  0.7356 0.9615 0.9898 1.00000

 # Gráfico de sedimentación (scree plot) para el análisis sin estandarización
ggscreeplot(nad.pca.cov, type = c('pev', 'cev')) +
  xlab('Número de componentes principales') +
  ylab('Proporción de la varianza explicada') +
  geom_line(colour = 'royalblue') +
  geom_point(colour = 'royalblue')

# Detalles adicionales sobre las componentes principales
print("Varianza explicada por cada componente (con estandarización):")

[1] "Varianza explicada por cada componente (con estandarización):"

print(summary(nad.pca.cor)$importance)

                            PC1       PC2       PC3       PC4
Standard deviation     1.709821 0.9573003 0.3483069 0.1969033
Proportion of Variance 0.730870 0.2291100 0.0303300 0.0096900
Cumulative Proportion  0.730870 0.9599800 0.9903100 1.0000000

# Gráfico de sedimentación (scree plot) para el PCA con estandarización
ggscreeplot(nad.pca.cor, type = c('pev', 'cev')) +
  xlab('Número de componentes principales') +
  ylab('Proporción de la varianza explicada') +
  geom_line(colour = 'royalblue') +
  geom_point(colour = 'royalblue') +
  ggtitle("Gráfico de sedimentación (Scree Plot)")

# Visualización del biplot usando ggbiplot
ggbiplot(nad.pca.cor, labels = rownames(nadadores), ellipse = TRUE, circle = TRUE) +
  ggtitle("Biplot de Componentes Principales (Estandarizado)") +
  theme_minimal()

library(factoextra)

Warning: package 'factoextra' was built under R version 4.4.3

Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

# Usar factoextra para gráficos avanzados
# Varianza acumulada explicada por cada componente
fviz_eig(nad.pca.cor, addlabels = TRUE, ylim = c(0, 60))+
  ggtitle("Varianza Explicada por Cada Componente")

library(factoextra)
# Biplot con factoextra
fviz_pca_biplot(nad.pca.cor, 
                col.ind = "cos2", # Color de individuos basado en calidad de representación
                gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
                col.var = "contrib", # Color de variables basado en contribución
                repel = TRUE) + # Evitar superposición de etiquetas
  ggtitle("Biplot Avanzado con Factoextra")

# Contribución de las variables a las componentes principales
fviz_contrib(nad.pca.cor, choice = "var", axes = 1) +
  ggtitle("Contribución de Variables a la PC1")

# Contribución de las variables a las componentes principales
fviz_contrib(nad.pca.cor, choice = "var", axes = 2) +
  ggtitle("Contribución de Variables a la PC2")

# Calidad de representación de los individuos en el espacio de las componentes
fviz_pca_ind(nad.pca.cor, 
             col.ind = "cos2", # Color basado en calidad de representación
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE) +
  ggtitle("Calidad de Representación de Individuos")

library(ggplot2)
library(gridExtra)
# Combinar múltiples gráficos en una sola figura
grid.arrange(
  fviz_eig(nad.pca.cor, addlabels = TRUE),
  fviz_pca_biplot(nad.pca.cor, col.ind = "cos2", col.var = "contrib", repel = TRUE),
  fviz_contrib(nad.pca.cor, choice = "var", axes = 1),
  ncol = 2
)

# Extraer la calidad de representación (cos²) de los individuos
cos2_pca_cor <- nad.pca.cor$ind$cos2 # Para PCA estandarizado
cos2_pca_cov <- nad.pca.cov$ind$cos2 # Para PCA no estandarizado

PCA CON FACTOMINER

# Instalar y cargar el paquete FactoMineR
#if (!requireNamespace("FactoMineR", quietly = TRUE)) {
 # install.packages("FactoMineR")}

library(FactoMineR)

Warning: package 'FactoMineR' was built under R version 4.4.3

# Realizar el PCA con FactoMineR
pca_factomine <- PCA(nadadores, scale.unit = TRUE, graph = FALSE)

# Resumen del PCA
summary(pca_factomine)


Call:
PCA(X = nadadores, scale.unit = TRUE, graph = FALSE) 


Eigenvalues
                       Dim.1   Dim.2   Dim.3   Dim.4
Variance               2.923   0.916   0.121   0.039
% of var.             73.087  22.911   3.033   0.969
Cumulative % of var.  73.087  95.998  99.031 100.000

Individuals (the 10 first)
        Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr
1   |  1.188 | -0.044  0.005  0.001 |  1.153 10.355  0.942 |  0.280  4.608
2   |  2.304 |  1.926  9.063  0.699 |  1.208 11.383  0.275 | -0.316  5.896
3   |  1.678 |  0.705  1.213  0.176 |  1.464 16.709  0.761 |  0.236  3.279
4   |  1.549 | -1.305  4.163  0.710 |  0.822  5.262  0.281 | -0.145  1.246
5   |  3.055 | -3.050 22.730  0.997 |  0.008  0.001  0.000 |  0.091  0.490
6   |  2.332 | -2.296 12.876  0.969 |  0.269  0.564  0.013 |  0.229  3.095
7   |  2.386 | -2.093 10.706  0.769 | -0.983  7.533  0.170 | -0.587 20.269
8   |  1.461 | -0.826  1.666  0.319 | -1.156 10.420  0.626 |  0.267  4.184
9   |  1.611 |  0.897  1.964  0.310 | -1.291 12.987  0.642 | -0.180  1.916
10  |  0.857 |  0.188  0.086  0.048 | -0.562  2.465  0.430 |  0.595 20.823
      cos2  
1    0.055 |
2    0.019 |
3    0.020 |
4    0.009 |
5    0.001 |
6    0.010 |
7    0.060 |
8    0.033 |
9    0.013 |
10   0.481 |

Variables
       Dim.1    ctr   cos2    Dim.2    ctr   cos2    Dim.3    ctr   cos2  
tr1 |  0.887 26.934  0.787 | -0.435 20.685  0.190 | -0.063  3.295  0.004 |
tr2 |  0.851 24.744  0.723 | -0.505 27.836  0.255 |  0.068  3.795  0.005 |
tr3 |  0.833 23.759  0.695 |  0.498 27.024  0.248 |  0.238 46.854  0.057 |
tr4 |  0.847 24.563  0.718 |  0.473 24.455  0.224 | -0.236 46.056  0.056 |

# Gráficos automáticos generados por FactoMineR
plot(pca_factomine, choix = "ind") # Gráfico de individuos

plot(pca_factomine, choix = "var") # Gráfico de variables

# Visualización avanzada con factoextra
library(factoextra)
fviz_pca_ind(pca_factomine, col.ind = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE) +
  ggtitle("Individuos en PCA (FactoMineR)")

fviz_pca_var(pca_factomine, col.var = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE) +
  ggtitle("Variables en PCA (FactoMineR)")

library(ade4)

Warning: package 'ade4' was built under R version 4.4.3


Adjuntando el paquete: 'ade4'

The following object is masked from 'package:FactoMineR':

    reconst

# Realizar el PCA con ade4
pca_ade4 <- dudi.pca(nadadores, scale = TRUE, scannf = FALSE, nf = 2)

# Resumen del PCA
pca_ade4

Duality diagramm
class: pca dudi
$call: dudi.pca(df = nadadores, scale = TRUE, scannf = FALSE, nf = 2)

$nf: 2 axis-components saved
$rank: 4
eigen values: 2.923 0.9164 0.1213 0.03877
  vector length mode    content       
1 $cw    4      numeric column weights
2 $lw    14     numeric row weights   
3 $eig   4      numeric eigen values  

  data.frame nrow ncol content             
1 $tab       14   4    modified array      
2 $li        14   2    row coordinates     
3 $l1        14   2    row normed scores   
4 $co        4    2    column coordinates  
5 $c1        4    2    column normed scores
other elements: cent norm

# Visualización avanzada con factoextra
library(factoextra)
fviz_pca_ind(pca_ade4, col.ind = "cos2", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE) +
  ggtitle("Individuos en PCA (ade4)")

fviz_pca_var(pca_ade4, col.var = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE) +
  ggtitle("Variables en PCA (ade4)")

PCA CON EL PAQUETE ade4

library(ade4)
acp <- dudi.pca(nadadores,
                scannf=FALSE, 
                scale=TRUE, #Matriz de correlación 
                nf=2)

summary(acp)

Class: pca dudi
Call: dudi.pca(df = nadadores, scale = TRUE, scannf = FALSE, nf = 2)

Total inertia: 4

Eigenvalues:
    Ax1     Ax2     Ax3     Ax4 
2.92349 0.91642 0.12132 0.03877 

Projected inertia (%):
    Ax1     Ax2     Ax3     Ax4 
73.0872 22.9106  3.0329  0.9693 

Cumulative projected inertia (%):
    Ax1   Ax1:2   Ax1:3   Ax1:4 
  73.09   96.00   99.03  100.00

print(acp)

Duality diagramm
class: pca dudi
$call: dudi.pca(df = nadadores, scale = TRUE, scannf = FALSE, nf = 2)

$nf: 2 axis-components saved
$rank: 4
eigen values: 2.923 0.9164 0.1213 0.03877
  vector length mode    content       
1 $cw    4      numeric column weights
2 $lw    14     numeric row weights   
3 $eig   4      numeric eigen values  

  data.frame nrow ncol content             
1 $tab       14   4    modified array      
2 $li        14   2    row coordinates     
3 $l1        14   2    row normed scores   
4 $co        4    2    column coordinates  
5 $c1        4    2    column normed scores
other elements: cent norm

# Valores propios
acp$eig

[1] 2.92348745 0.91642395 0.12131769 0.03877091

inertia.dudi(acp)

Inertia information:
Call: inertia.dudi(x = acp)

Decomposition of total inertia:
    inertia     cum  cum(%)
Ax1 2.92349   2.923   73.09
Ax2 0.91642   3.840   96.00
Ax3 0.12132   3.961   99.03
Ax4 0.03877   4.000  100.00

# Vectores propios
acp$c1

           CS1        CS2
tr1 -0.5189780  0.4548106
tr2 -0.4974341  0.5275946
tr3 -0.4874313 -0.5198444
tr4 -0.4956126 -0.4945231

# Primera forma 
plot(acp$eig,type="b",pch=20,col="blue")
abline(h=1,lty=3,col="red")

# Segunda forma
screeplot(acp, main ="Screeplot - Valores Propios")

# Tercera forma
library(factoextra)
eig.val <- get_eigenvalue(acp)
eig.val

      eigenvalue variance.percent cumulative.variance.percent
Dim.1 2.92348745       73.0871862                    73.08719
Dim.2 0.91642395       22.9105988                    95.99779
Dim.3 0.12131769        3.0329421                    99.03073
Dim.4 0.03877091        0.9692728                   100.00000

barplot(eig.val[, 2], names.arg=1:nrow(eig.val), 
        main = "Autovalores",
        xlab = "Componentes Principales",
        ylab = "Porcentaje de variancias",
        col ="steelblue")
lines(x = 1:nrow(eig.val), eig.val[, 2], 
      type="b", pch=19, col = "red")

# Cuarta forma
library(factoextra)
fviz_screeplot(acp, ncp=6)

fviz_eig(acp, addlabels=TRUE, hjust = -0.3)

fviz_eig(acp, addlabels=TRUE, hjust = -0.3,
         barfill="white", barcolor ="darkblue",
         linecolor ="red") + ylim(0,50) + theme_minimal()

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_bar()`).

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_line()`).

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_text()`).

# Correlaciones entre las variables y los componentes
acp$co

         Comp1      Comp2
tr1 -0.8873595  0.4353904
tr2 -0.8505232  0.5050665
tr3 -0.8334202 -0.4976472
tr4 -0.8474088 -0.4734071

# Gráfica de Variables sobre el círculo de correlaciones
# Primera forma 
s.corcircle(acp$co,grid=FALSE)

# Segunda forma
fviz_pca_var(acp)

fviz_pca_var(acp, col.var="steelblue")+theme_minimal()

#Contribución de las variables a los componentes principaes

(contrib <- acp$co*acp$co)

        Comp1     Comp2
tr1 0.7874068 0.1895648
tr2 0.7233898 0.2550922
tr3 0.6945892 0.2476527
tr4 0.7181017 0.2241143

contrib <- as.matrix(contrib)

library(corrplot)

Warning: package 'corrplot' was built under R version 4.4.2

corrplot 0.95 loaded

corrplot(contrib,is.corr=FALSE)

# Scores o Puntuaciones de cada individuo
acp$li[1:10,]

         Axis1        Axis2
1   0.04400452 -1.152644418
2  -1.92593589 -1.208460292
3  -0.70459158 -1.464136885
4   1.30535238 -0.821677830
5   3.05013379 -0.008340468
6   2.29566526 -0.268999975
7   2.09327497  0.983082891
8   0.82571616  1.156261799
9  -0.89663506  1.290798352
10 -0.18770397  0.562335102

#Análisis descriptivo de los scores

options(scipen=999)
cov(acp$li)

                        Axis1                   Axis2
Axis1 3.148371098924580557821 0.000000000000001906836
Axis2 0.000000000000001906836 0.986918103435560745140

cor(acp$li)

                        Axis1                   Axis2
Axis1 1.000000000000000000000 0.000000000000001081758
Axis2 0.000000000000001081758 1.000000000000000000000

library(psych)

Warning: package 'psych' was built under R version 4.4.3


Adjuntando el paquete: 'psych'

The following objects are masked from 'package:scales':

    alpha, rescale

The following objects are masked from 'package:ggplot2':

    %+%, alpha

describe(acp$li)

      vars  n mean   sd median trimmed  mad   min  max range  skew kurtosis
Axis1    1 14    0 1.77  -0.07   -0.03 1.69 -2.67 3.05  5.72  0.08    -1.19
Axis2    2 14    0 0.99   0.23    0.01 1.31 -1.46 1.29  2.75 -0.15    -1.71
        se
Axis1 0.47
Axis2 0.27

# Gráfica de individuos sobre el primer plano de componentes
# Primera forma
s.label(acp$li,xax=1,yax=2,clabel=0.7,grid=FALSE,boxes=FALSE)

# Segunda forma 
fviz_pca_ind(acp)

# Gráfica de individuos sobre los componentes 2 y 3. Para realizar este gráfico se tiene que cambiar el número de factores en la función dudi.pca
# s.label(acp$li,xax=2,yax=3,clabel=0.7,grid=FALSE,boxes=FALSE)

# Gráfica de individuos sobre el primer plano con biplot
# Primera forma
s.label(acp$li,clabel=0.7,grid=FALSE,boxes=FALSE)
s.corcircle(acp$co,grid=FALSE,add=TRUE,clabel=0.7)

# Segunda forma
fviz_pca_biplot(acp, repel = FALSE,
                col.var = "steelblue",
                col.ind = "black" )

#Grabar los datos y los resultados de los scores en un archivo CSV

salida.acp <- cbind(nadadores,acp$li[,c(1,2)])
head(salida.acp)

  tr1 tr2 tr3 tr4       Axis1        Axis2
1  10  10  13  12  0.04400452 -1.152644418
2  12  12  14  15 -1.92593589 -1.208460292
3  11  10  14  13 -0.70459158 -1.464136885
4   9   9  11  11  1.30535238 -0.821677830
5   8   8   9   8  3.05013379 -0.008340468
6   8   9  10   9  2.29566526 -0.268999975

write.csv(salida.acp,"nad-acp.csv")

import pandas as pd
ruta_archivo = "nadadores.xlsx"

# Cargar el archivo
nadadores = pd.read_excel(ruta_archivo)

# Mostrar el DataFrame original
print("DataFrame original:")

DataFrame original:

print(nadadores.head())

   Nadador  tr1  tr2  tr3  tr4
0        1   10   10   13   12
1        2   12   12   14   15
2        3   11   10   14   13
3        4    9    9   11   11
4        5    8    8    9    8

# Eliminar la primera columna (por posición)
nadadores = nadadores.iloc[:, 1:]

# Mostrar el DataFrame modificado
print("\nDataFrame después de eliminar la primera columna:")


DataFrame después de eliminar la primera columna:

print(nadadores.head())

   tr1  tr2  tr3  tr4
0   10   10   13   12
1   12   12   14   15
2   11   10   14   13
3    9    9   11   11
4    8    8    9    8

PCA EN PYTHON

# 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')

# Resumen de la base de datos
print(nadadores.info())  # Tipos de datos y valores no nulos

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 14 entries, 0 to 13
Data columns (total 4 columns):
 #   Column  Non-Null Count  Dtype
---  ------  --------------  -----
 0   tr1     14 non-null     int64
 1   tr2     14 non-null     int64
 2   tr3     14 non-null     int64
 3   tr4     14 non-null     int64
dtypes: int64(4)
memory usage: 580.0 bytes
None

print(nadadores.describe())  # Estadísticas descriptivas

             tr1        tr2        tr3        tr4
count  14.000000  14.000000  14.000000  14.000000
mean   11.214286  11.214286  11.571429  11.285714
std     2.224983   2.082106   1.910066   2.127786
min     8.000000   8.000000   8.000000   8.000000
25%    10.000000  10.000000  10.250000   9.250000
50%    11.000000  11.000000  11.500000  11.000000
75%    12.750000  12.750000  13.000000  13.000000
max    15.000000  15.000000  14.000000  15.000000

print('----------------------')

----------------------

print('Media de cada variable')

Media de cada variable

print('----------------------')

----------------------

nadadores.mean(axis=0)

tr1    11.214286
tr2    11.214286
tr3    11.571429
tr4    11.285714
dtype: float64

print('-------------------------')

-------------------------

print('Varianza de cada variable')

Varianza de cada variable

print('-------------------------')

-------------------------

nadadores.var(axis=0)

tr1    4.950549
tr2    4.335165
tr3    3.648352
tr4    4.527473
dtype: float64

# Entrenamiento modelo PCA con escalado de los datos
# ==============================================================================
pca_pipe = make_pipeline(StandardScaler(), PCA())
pca_pipe.fit(nadadores)

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

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


# Se extrae el modelo entrenado del pipeline
modelo_pca = pca_pipe.named_steps['pca']

# Se convierte el array a dataframe para añadir nombres a los ejes.
pd.DataFrame(
    data    = modelo_pca.components_,
    columns = nadadores.columns,
    index   = ['PC1', 'PC2', 'PC3', 'PC4']
)

          tr1       tr2       tr3       tr4
PC1  0.518978  0.497434  0.487431  0.495613
PC2  0.454811  0.527595 -0.519844 -0.494523
PC3 -0.181519  0.194812  0.684497 -0.678648
PC4  0.700614 -0.660494  0.153742 -0.221927

#La influencia de las variables en cada componente 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(nadadores.columns)), nadadores.columns)

([<matplotlib.axis.YTick object at 0x0000023C3E1137A0>, <matplotlib.axis.YTick object at 0x0000023C3E137020>, <matplotlib.axis.YTick object at 0x0000023C3E0E3110>, <matplotlib.axis.YTick object at 0x0000023C3E152E40>], [Text(0, 0, 'tr1'), Text(0, 1, 'tr2'), Text(0, 2, 'tr3'), Text(0, 3, 'tr4')])

plt.xticks(range(len(nadadores.columns)), np.arange(modelo_pca.n_components_) + 1)

([<matplotlib.axis.XTick object at 0x0000023C3E134B30>, <matplotlib.axis.XTick object at 0x0000023C3E136FF0>, <matplotlib.axis.XTick object at 0x0000023C3E153B60>, <matplotlib.axis.XTick object at 0x0000023C3E198620>], [Text(0, 0, '1'), Text(1, 0, '2'), Text(2, 0, '3'), Text(3, 0, '4')])

plt.grid(False)
plt.colorbar();
plt.show()

# Porcentaje de varianza explicada por cada componente
# ==============================================================================
print('----------------------------------------------------')

----------------------------------------------------

print('Porcentaje de varianza explicada por cada componente')

Porcentaje de varianza explicada por cada componente

print('----------------------------------------------------')

----------------------------------------------------

print(modelo_pca.explained_variance_ratio_)

[0.73087186 0.22910599 0.03032942 0.00969273]

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(6, 4))
ax.bar(
    x      = np.arange(modelo_pca.n_components_) + 1,
    height = modelo_pca.explained_variance_ratio_
)

for x, y in zip(np.arange(len(nadadores.columns)) + 1, modelo_pca.explained_variance_ratio_):
    label = round(y, 2)
    ax.annotate(
        label,
        (x,y),
        textcoords="offset points",
        xytext=(0,10),
        ha='center'
    )

ax.set_xticks(np.arange(modelo_pca.n_components_) + 1)
ax.set_ylim(0, 1.1)

(0.0, 1.1)

ax.set_title('Porcentaje de varianza explicada por cada componente')
ax.set_xlabel('Componente principal')
ax.set_ylabel('Por. varianza explicada');
plt.show()

# Porcentaje de varianza explicada acumulada
# ==============================================================================
prop_varianza_acum = modelo_pca.explained_variance_ratio_.cumsum()
print('------------------------------------------')

------------------------------------------

print('Porcentaje de varianza explicada acumulada')

Porcentaje de varianza explicada acumulada

print('------------------------------------------')

------------------------------------------

print(prop_varianza_acum)

[0.73087186 0.95997785 0.99030727 1.        ]

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(6, 4))
ax.plot(
    np.arange(len(nadadores.columns)) + 1,
    prop_varianza_acum,
    marker = 'o'
)

for x, y in zip(np.arange(len(nadadores.columns)) + 1, prop_varianza_acum):
    label = round(y, 2)
    ax.annotate(
        label,
        (x,y),
        textcoords="offset points",
        xytext=(0,10),
        ha='center'
    )
    
ax.set_ylim(0, 1.1)

(0.0, 1.1)

ax.set_xticks(np.arange(modelo_pca.n_components_) + 1)
ax.set_title('Porcentaje de varianza explicada acumulada')
ax.set_xlabel('Componente principal')
ax.set_ylabel('Por. varianza acumulada');
plt.show()

#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=nadadores)
proyecciones = pd.DataFrame(
    proyecciones,
    columns = ['PC1', 'PC2', 'PC3', 'PC4'],
    index   = nadadores.index
)
proyecciones.head()

        PC1       PC2       PC3       PC4
0 -0.044005 -1.152644  0.279755  0.044962
1  1.925936 -1.208460 -0.316442 -0.201074
2  0.704592 -1.464137  0.235998  0.347026
3 -1.305352 -0.821678 -0.145476 -0.011432
4 -3.050134 -0.008340  0.091265  0.148648

#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(nadadores).T)
proyecciones = pd.DataFrame(proyecciones, index = ['PC1', 'PC2', 'PC3', 'PC4'])
proyecciones = proyecciones.transpose().set_index(nadadores.index)
proyecciones.head()

        PC1       PC2       PC3       PC4
0 -0.044005 -1.152644  0.279755  0.044962
1  1.925936 -1.208460 -0.316442 -0.201074
2  0.704592 -1.464137  0.235998  0.347026
3 -1.305352 -0.821678 -0.145476 -0.011432
4 -3.050134 -0.008340  0.091265  0.148648

#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

# Recostruccion de las proyecciones
# ==============================================================================
recostruccion = pca_pipe.inverse_transform(X=proyecciones)
recostruccion = pd.DataFrame(
                    recostruccion,
                    columns = nadadores.columns,
                    index   = nadadores.index
)
print('------------------')

------------------

print('Valores originales')

Valores originales

print('------------------')

------------------

print(recostruccion.head())

    tr1   tr2   tr3   tr4
0  10.0  10.0  13.0  12.0
1  12.0  12.0  14.0  15.0
2  11.0  10.0  14.0  13.0
3   9.0   9.0  11.0  11.0
4   8.0   8.0   9.0   8.0

print('---------------------')

---------------------

print('Valores reconstruidos')

Valores reconstruidos

print('---------------------')

---------------------

print(nadadores.head())

   tr1  tr2  tr3  tr4
0   10   10   13   12
1   12   12   14   15
2   11   10   14   13
3    9    9   11   11
4    8    8    9    8

# Verificar valores faltantes
print(nadadores.isnull().sum())

tr1    0
tr2    0
tr3    0
tr4    0
dtype: int64


# Imputar valores faltantes (ejemplo: media)
#nadadores.fillna(nadadores.mean(), inplace=True)

from sklearn.preprocessing import StandardScaler

# Seleccionar solo variables numéricas
variables_numericas = nadadores.select_dtypes(include=['float64', 'int64'])

# Estandarizar los datos
scaler = StandardScaler()
datos_estandarizados = scaler.fit_transform(variables_numericas)

from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

# Aplicar PCA
pca = PCA(n_components=2)  # Reducir a 2 componentes principales
componentes_principales = pca.fit_transform(datos_estandarizados)

# Convertir a DataFrame para mejor visualización
pca_df = pd.DataFrame(componentes_principales, columns=["PC1", "PC2"])
print(pca_df.head())

        PC1       PC2
0 -0.044005 -1.152644
1  1.925936 -1.208460
2  0.704592 -1.464137
3 -1.305352 -0.821678
4 -3.050134 -0.008340

# Varianza explicada
varianza_explicada = pca.explained_variance_ratio_
print("Varianza explicada por cada componente:", varianza_explicada)

Varianza explicada por cada componente: [0.73087186 0.22910599]

# Gráfico de dispersión de los datos en PC1 y PC2
plt.figure(figsize=(8, 6))
plt.scatter(pca_df["PC1"], pca_df["PC2"], c='blue', edgecolor='k', s=100)
plt.xlabel("Componente Principal 1")
plt.ylabel("Componente Principal 2")
plt.title("PCA - Proyección en las dos primeras componentes")
plt.grid(True)
plt.show()

# Convertir a DataFrame para mejor visualización
scores_df = pd.DataFrame(componentes_principales, columns=["PC1", "PC2"])
print("Scores (coordenadas de los individuos en las componentes principales):")

Scores (coordenadas de los individuos en las componentes principales):

print(scores_df.head())

        PC1       PC2
0 -0.044005 -1.152644
1  1.925936 -1.208460
2  0.704592 -1.464137
3 -1.305352 -0.821678
4 -3.050134 -0.008340

# Loadings
loadings = pd.DataFrame(pca.components_, columns=variables_numericas.columns, index=["PC1", "PC2"])
print("Loadings:")

Loadings:

print(loadings)

          tr1       tr2       tr3       tr4
PC1  0.518978  0.497434  0.487431  0.495613
PC2  0.454811  0.527595 -0.519844 -0.494523

# Biplot
plt.figure(figsize=(10, 8))
plt.scatter(scores_df["PC1"], scores_df["PC2"], color="blue", alpha=0.7, label="Individuos")

for i, feature in enumerate(variables_numericas.columns):
    plt.arrow(0, 0, pca.components_[0, i], pca.components_[1, i], color='red', head_width=0.05)
    plt.text(pca.components_[0, i] * 1.1, pca.components_[1, i] * 1.1, feature, color='red')

plt.axhline(0, color='gray', linestyle='--', linewidth=0.8)
plt.axvline(0, color='gray', linestyle='--', linewidth=0.8)
plt.xlabel("Componente Principal 1")
plt.ylabel("Componente Principal 2")
plt.title("Biplot: Scores y Loadings")
plt.legend()
plt.grid(True)
plt.show()

from sklearn.decomposition import PCA
import pandas as pd

# Supongamos que ya tenemos los datos estandarizados en 'datos_estandarizados'
pca = PCA(n_components=2)
pca.fit(datos_estandarizados)

PCA(n_components=2)

# Obtener los loadings
loadings = pd.DataFrame(pca.components_, columns=nadadores.columns, index=["PC1", "PC2"])
print("Loadings (contribución de cada variable a las componentes principales):")

Loadings (contribución de cada variable a las componentes principales):

print(loadings)

          tr1       tr2       tr3       tr4
PC1  0.518978  0.497434  0.487431  0.495613
PC2  0.454811  0.527595 -0.519844 -0.494523