A PCA (Principal Components Analysis ou Análise de Componentes Principais) é uma técnica de aprendizado não supervisionado amplamente utilizada para exploração de dados e redução de dimensionalidade. Seu principal objetivo é identificar a representação n-dimensional mais enxuta que preserve a maior variabilidade dos dados, especialmente em conjuntos com um grande número de variáveis (atributos, parâmetros ou dimensões). Além disso, a PCA também pode ser aplicada em clustering, que consiste em agrupar dados com base em padrões similares, facilitando sua categorização em conjuntos distintos.
Neste artigo, demonstramos como desenvolver a PCA manualmente no R, desde a preparação dos dados até a visualização dos resultados. Calculamos simultaneamente os eigenvalues ou autovalores e os eigenvectors ou autovetores de uma matriz quadrada, por meio da função eigen(). O prefixo alemão eigen pode ser traduzido aproximadamente como próprio ou característico. Um autovetor de A é um vetor que é transformado em um múltiplo de si mesmo pela transformação matricial T(x) = Ax. Por outro lado, podemos pensar em um autovetor como algo que descreve uma propriedade intrínseca, ou característica, de A.
Em álgebra linear, eigenvectors são vetores que mantêm sua direção após uma transformação linear, enquanto eigenvalues representam o fator pelo qual esses vetores são escalados. Na PCA, os eigenvectors representam os componentes principais do conjunto de dados, enquanto seus eigenvalues correspondentes indicam a quantidade de variância explicada por cada componente.
O processo manual da PCA envolve verificação de correlação (Bartlett), padronização dos dados, decomposição da matriz de covariância, cálculo de autovalores/vetores e projeção dos dados nos componentes.
Definimos o diretório de trabalho:
setwd("/home/gf/Scripts/Moyses")
Carregamos as bibliotecas necessárias:
O dataset contém 210 observações de 7 variáveis numéricas, que representam as medições das propriedades geométricas de grãos pertencentes a três variedades diferentes de trigo.
dados <- read.table('seeds_dataset.txt', header=TRUE, sep=' ')
Verificamos o tipo dos dados, se há registros faltantes e fornecemos a estatística básica.
str(dados)
## 'data.frame': 210 obs. of 7 variables:
## $ area : num 15.3 14.9 14.3 13.8 16.1 ...
## $ perimeter : num 14.8 14.6 14.1 13.9 15 ...
## $ compactness : num 0.871 0.881 0.905 0.895 0.903 ...
## $ length : num 5.76 5.55 5.29 5.32 5.66 ...
## $ width : num 3.31 3.33 3.34 3.38 3.56 ...
## $ asymmetry : num 2.22 1.02 2.7 2.26 1.35 ...
## $ length_grove: num 5.22 4.96 4.83 4.8 5.17 ...
colSums(is.na(dados))
## area perimeter compactness length width asymmetry
## 0 0 0 0 0 0
## length_grove
## 0
head(dados)
## area perimeter compactness length width asymmetry length_grove
## 1 15.26 14.84 0.8710 5.763 3.312 2.221 5.220
## 2 14.88 14.57 0.8811 5.554 3.333 1.018 4.956
## 3 14.29 14.09 0.9050 5.291 3.337 2.699 4.825
## 4 13.84 13.94 0.8955 5.324 3.379 2.259 4.805
## 5 16.14 14.99 0.9034 5.658 3.562 1.355 5.175
## 6 14.38 14.21 0.8951 5.386 3.312 2.462 4.956
summary(dados)
## area perimeter compactness length
## Min. :10.59 Min. :12.41 Min. :0.8081 Min. :4.899
## 1st Qu.:12.27 1st Qu.:13.45 1st Qu.:0.8569 1st Qu.:5.262
## Median :14.36 Median :14.32 Median :0.8734 Median :5.524
## Mean :14.85 Mean :14.56 Mean :0.8710 Mean :5.629
## 3rd Qu.:17.30 3rd Qu.:15.71 3rd Qu.:0.8878 3rd Qu.:5.980
## Max. :21.18 Max. :17.25 Max. :0.9183 Max. :6.675
## width asymmetry length_grove
## Min. :2.630 Min. :0.7651 Min. :4.519
## 1st Qu.:2.944 1st Qu.:2.5615 1st Qu.:5.045
## Median :3.237 Median :3.5990 Median :5.223
## Mean :3.259 Mean :3.7002 Mean :5.408
## 3rd Qu.:3.562 3rd Qu.:4.7687 3rd Qu.:5.877
## Max. :4.033 Max. :8.4560 Max. :6.550
O objetivo do Teste de Bartlett, também chamado Teste de Esfericidade de Bartlett, é assegurar que a matriz de correlação das variáveis seja significativamente divergente da Matriz Identidade, de modo que uma técnica de redução de dimensionalidade seja aplicável. Se o valor de p-value for menor que o nível de significância adotado (0.10, 0.05 ou 0.01), então o conjunto de dados avaliado pode ser submetido a alguma técnica de redução de dimensionalidade.
bart <- function(dat){
R <- cor(dados)
p <- ncol(dados)
n <- nrow(dados)
chi2 <- -((n - 1)-((2 * p)+ 5)/ 6) * log(det(R))
df <- (p * (p - 1)/ 2)
crit <- qchisq(.95, df) # Nível de significância
p <- pchisq(chi2, df, lower.tail = F) # p-value
cat("Bartlett's test: X2(",
df,") = ", chi2,", p = ",
round(p, 5),
sep="" )
}
bart(dados)
## Bartlett's test: X2(21) = 3623.407, p = 0
Nesse caso, o p-valor é < .001, indicando que as variáveis são altamente correlacionadas e, portanto, a PCA é apropriada.
Como as variáveis têm escalas diferentes, padronizamos os dados para média 0 e desvio padrão 1 usando a função scale(). Isso evita que variáveis com maior variância dominem a análise.
dados_norm <- scale(dados)
Representação matricial das observações
Transformamos a tabela em uma matriz de dados:
m1 <- as.matrix(dados_norm)
head(m1)
## area perimeter compactness length width asymmetry
## 1 0.14175904 0.214948819 6.045733e-05 0.30349301 0.1413640 -0.9838010
## 2 0.01116136 0.008204153 4.274938e-01 -0.16822270 0.1969616 -1.7839036
## 3 -0.19160873 -0.359341919 1.438945e+00 -0.76181710 0.2075516 -0.6658882
## 4 -0.34626388 -0.474200066 1.036904e+00 -0.68733567 0.3187467 -0.9585276
## 5 0.44419577 0.329806966 1.371233e+00 0.06650665 0.8032397 -1.5597684
## 6 -0.16067770 -0.267455401 1.019976e+00 -0.54740087 0.1413640 -0.8235144
## length_grove
## 1 -0.3826631
## 2 -0.9198156
## 3 -1.1863572
## 4 -1.2270506
## 5 -0.4742231
## 6 -0.9198156
Calculamos o vetor de médias amostral
media_m1 <- round(apply(m1, 2, mean))
media_m1
## area perimeter compactness length width asymmetry
## 0 0 0 0 0 0
## length_grove
## 0
Calculamos o vetor de desvio padrão amostral
dp_m1 <- apply(m1, 2, sd)
dp_m1
## area perimeter compactness length width asymmetry
## 1 1 1 1 1 1
## length_grove
## 1
Calculamos a matriz de covariâncias amostral
A matriz de covariância é uma matriz quadrada que resume as covariâncias entre um conjunto de variáveis.
cov_m1 <- cov(m1)
cov_m1
## area perimeter compactness length width
## area 1.0000000 0.9943409 0.6082884 0.9499854 0.9707706
## perimeter 0.9943409 1.0000000 0.5292436 0.9724223 0.9448294
## compactness 0.6082884 0.5292436 1.0000000 0.3679151 0.7616345
## length 0.9499854 0.9724223 0.3679151 1.0000000 0.8604149
## width 0.9707706 0.9448294 0.7616345 0.8604149 1.0000000
## asymmetry -0.2295723 -0.2173404 -0.3314709 -0.1715624 -0.2580365
## length_grove 0.8636927 0.8907839 0.2268248 0.9328061 0.7491315
## asymmetry length_grove
## area -0.22957233 0.86369275
## perimeter -0.21734037 0.89078390
## compactness -0.33147087 0.22682482
## length -0.17156243 0.93280609
## width -0.25803655 0.74913147
## asymmetry 1.00000000 -0.01107902
## length_grove -0.01107902 1.00000000
Identificamos as variáveis correlacionadas
Executamos a função corr_check(), que imprime automaticamente os pares de variáveis altamente correlacionadas.
corr_check <- function(Dataset, threshold){
matriz_cor <- cor(Dataset)
matriz_cor
for (i in 1:nrow(matriz_cor)){
correlations <- which((abs(matriz_cor[i,i:ncol(matriz_cor)]) > threshold) & (matriz_cor[i,i:ncol(matriz_cor)] != 1))
if(length(correlations)> 0){
lapply(correlations,FUN = function(x) (cat(paste(colnames(Dataset)[i], "com",colnames(Dataset)[x]), "\n")))
}
}
}
corr_check(dados, 0.85)
## area com perimeter
## area com length
## area com width
## area com length_grove
## perimeter com compactness
## perimeter com length
## perimeter com asymmetry
## length com perimeter
## length com length
A decomposição espectral da matriz retorna uma lista com dois componentes: um vetor contendo os eigenvalues e uma matriz cujas colunas contêm os engenvectors da matriz de covariância.
autos <- eigen(cov_m1)
autos
## eigen() decomposition
## $values
## [1] 5.0312011860 1.1975728470 0.6780034386 0.0683644770 0.0187136090
## [6] 0.0053320457 0.0008123968
##
## $vectors
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -0.4444735 0.02656355 -0.02587094 0.19363997 -0.20441167 -0.42643686
## [2,] -0.4415715 0.08400282 0.05983912 0.29545659 -0.17427591 -0.47623853
## [3,] -0.2770174 -0.52915125 -0.62969178 -0.33281640 0.33265481 -0.14162884
## [4,] -0.4235633 0.20597518 0.21187966 0.26340659 0.76609839 0.27357647
## [5,] -0.4328187 -0.11668963 -0.21648338 0.19963039 -0.46536555 0.70301171
## [6,] 0.1186925 0.71688203 -0.67950584 0.09246481 0.03625822 -0.01964186
## [7,] -0.3871608 0.37719327 0.21389720 -0.80414995 -0.11134657 0.04282974
## [,7]
## [1,] 0.734805689
## [2,] -0.670751532
## [3,] -0.072552703
## [4,] 0.046276051
## [5,] -0.039289079
## [6,] -0.003723456
## [7,] -0.034498098
Matriz de variâncias dos componentes principais
A matriz V é construída para armazenar os autovalores em sua diagonal, representando as variâncias dos componentes principais. Essa matriz é fundamental para reconstruir a matriz de covariância original e calcular os loadings ajustados, que representam as direções desses componentes no espaço original. Sua diagonal confirma que os componentes são ortogonais.
L <- autos$values # Armazena os autovalores em L
V <- diag(L) # Insere os autovalores na diagonal de V
V # Exibe a matriz com autovalores na diagonal
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 5.031201 0.000000 0.0000000 0.00000000 0.00000000 0.000000000 0.0000000000
## [2,] 0.000000 1.197573 0.0000000 0.00000000 0.00000000 0.000000000 0.0000000000
## [3,] 0.000000 0.000000 0.6780034 0.00000000 0.00000000 0.000000000 0.0000000000
## [4,] 0.000000 0.000000 0.0000000 0.06836448 0.00000000 0.000000000 0.0000000000
## [5,] 0.000000 0.000000 0.0000000 0.00000000 0.01871361 0.000000000 0.0000000000
## [6,] 0.000000 0.000000 0.0000000 0.00000000 0.00000000 0.005332046 0.0000000000
## [7,] 0.000000 0.000000 0.0000000 0.00000000 0.00000000 0.000000000 0.0008123968
Para reproduzir a matriz de correlação original, podemos utililizar os eigenvalues e eigenvectors.
autos$vectors %*% V %*% t(autos$vectors) # V L V
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1.0000000 0.9943409 0.6082884 0.9499854 0.9707706 -0.22957233
## [2,] 0.9943409 1.0000000 0.5292436 0.9724223 0.9448294 -0.21734037
## [3,] 0.6082884 0.5292436 1.0000000 0.3679151 0.7616345 -0.33147087
## [4,] 0.9499854 0.9724223 0.3679151 1.0000000 0.8604149 -0.17156243
## [5,] 0.9707706 0.9448294 0.7616345 0.8604149 1.0000000 -0.25803655
## [6,] -0.2295723 -0.2173404 -0.3314709 -0.1715624 -0.2580365 1.00000000
## [7,] 0.8636927 0.8907839 0.2268248 0.9328061 0.7491315 -0.01107902
## [,7]
## [1,] 0.86369275
## [2,] 0.89078390
## [3,] 0.22682482
## [4,] 0.93280609
## [5,] 0.74913147
## [6,] -0.01107902
## [7,] 1.00000000
Calculamos o percentual de explicação, através da divisão dos eigenvalues pela soma da diagonal da matriz de convariâncias:
Calculamos a importância dos Componentes
Calculamos a proporção de variância explicada por cada componente.
pexp <- autos$values/sum(diag(cov_m1))
pexp
## [1] 0.7187430266 0.1710818353 0.0968576341 0.0097663539 0.0026733727
## [6] 0.0007617208 0.0001160567
Calculamos a proporção de variância acumulada, através da função cumsum()
pexp_ac <- cumsum(pexp)
pexp_ac
## [1] 0.7187430 0.8898249 0.9866825 0.9964488 0.9991222 0.9998839 1.0000000
Embora o primeiro componente seja responsável pela explicação de mais de 70% da variância observada, escolhemos os três primeiros para descrever a variância total contida nos dados.
pca.escores <- m1 %*% autos$vectors
colnames(pca.escores) <- c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7")
head(pca.escores)
## PC1 PC2 PC3 PC4 PC5 PC6
## 1 -0.316291277 -0.7818009 0.62950581 0.41585226 0.107239413 0.02251669
## 2 0.003378106 -1.9086530 0.66815790 0.43203823 -0.044300562 0.01887606
## 3 0.458348123 -1.9026788 -0.93026596 0.11102631 0.008202305 -0.05119072
## 4 0.590524599 -1.9264659 -0.49812067 0.22843077 -0.074674108 0.22895337
## 5 -1.100280938 -2.0631601 -0.05657024 0.14208115 -0.018729812 0.05251314
## 6 0.336012252 -1.6330142 -0.43786687 0.09795543 0.006165489 -0.02216203
## PC7
## 1 -0.014661607
## 2 -0.005466297
## 3 -0.004167310
## 4 -0.010026910
## 5 -0.000620942
## 6 -0.008760160
Padronizamos os sinais dos autovetores
Um problema crítico no cálculo manual da PCA é a inversão arbitrária de sinais em autovetores, detectada quando confrontados com o resultado das funções prcomp()/princomp().
Para padronização consistente dos autovetores, garantimos que a soma dos elementos de cada autovetor seja positiva:
for (i in 1:ncol(autos$vectors)) {
if (sum(autos$vectors[, i]) < 0) {
autos$vectors[, i] <- autos$vectors[, i] * -1
pca.escores[, i] <- pca.escores[, i] * -1
}
}
pca = list("loadings" = autos$vectors, "scores" = pca.escores)
# Ajuste das margens (bottom, left, top, right) e tamanho da fonte
par(mar = c(1, 1, 1, 1) + 0.1, cex = 1.0)
# Plotar o biplot
biplot(pca$scores, pca$loadings,
xlim = c(-4, 5),
ylim = c(-3, 3),
cex = c(0.8, 1.2)) # cex: tamanho dos textos (escores e loadings)
Observamos, no gráfico, uma leve correlação entre as variáveis length grove, length, perimeter, área e width com a PC1 e uma acentuada correlação entre a variável asymmetry e a PC2. As variáveis correlacionadas com a PC1 estão positivamente correlacionadas entre si, enquanto as variáveis asymmetry e compactness possuem correlação negativa.
O gráfico de dispersão dos dois primeiros componentes ajuda a identificar padrões ou agrupamentos nos dados:
plot(x = pca.escores[, 1], y= pca.escores[, 2],
pch = c(rep(21, 70), rep(22, 70), rep(24, 70)),
col = c(rep("blue", "70"), rep("red", 70), rep("green", 70)),
xlab = "PC1", ylab = "PC2", main = "Primeiro x Segundo Componente Principal",
cex.lab = 1.2, cex.axis = 1.2)
round(cor(pca.escores), 2)
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## PC1 1 0 0 0 0 0 0
## PC2 0 1 0 0 0 0 0
## PC3 0 0 1 0 0 0 0
## PC4 0 0 0 1 0 0 0
## PC5 0 0 0 0 1 0 0
## PC6 0 0 0 0 0 1 0
## PC7 0 0 0 0 0 0 1
Os loadings ajustados são calculados multiplicando os autovetores pela raiz quadrada dos autovalores, representando as correlações entre variáveis originais e componentes principais. Isso permite interpretar a contribuição de cada variável para a redução de dimensionalidade.
loadings <- autos$vectors %*% sqrt(V)
colnames(loadings) <- paste0('PC', 1:7)
loadings
## PC1 PC2 PC3 PC4 PC5 PC6
## [1,] 0.9969692 0.02906947 0.02130238 -0.05063027 -0.02796304 0.031138785
## [2,] 0.9904598 0.09192737 -0.04927210 -0.07725186 -0.02384054 0.034775346
## [3,] 0.6213594 -0.57906965 0.51849428 0.08702018 0.04550641 0.010341859
## [4,] 0.9500669 0.22540620 -0.17446376 -0.06887187 0.10480049 -0.019976789
## [5,] 0.9708269 -0.12769775 0.17825451 -0.05219656 -0.06366093 -0.051334518
## [6,] -0.2662313 0.78451032 0.55951166 -0.02417641 0.00496004 0.001434266
## [7,] 0.8684149 0.41277645 -0.17612501 0.21025788 -0.01523195 -0.003127464
## PC7
## [1,] -0.0209438545
## [2,] 0.0191181460
## [3,] 0.0020679389
## [4,] -0.0013189866
## [5,] 0.0011198399
## [6,] 0.0001061281
## [7,] 0.0009832846
Considerações finais
Este artigo demonstrou como implementar a PCA passo a passo no R, desde a análise de correlação até a visualização dos componentes principais.
A função princomp() utiliza a abordagem computacional denominada decomposição espectral ou eigendecomposition. Para uma implementação mais eficiente, utilize a função prcomp(), que adota a abordagem chamada singular value decomposition (SVD) ou decomposição em valores singulares.
Referências
Boehmke, B. (2016). Principal Components Analysis. In UC Business Analytics R Programming Guide. University of Cincinnati. Disponível em: https://uc-r.github.io/pca
Charytanowicz, M., Niewczas, J., Kulczycki, P., Kowalski, P., & Lukasik, S. (2010). Seeds [Dataset]. UCI Machine Learning Repository. https://doi.org/10.24432/C5H30K.
Hartmann, K., Krois, J., Rudolph, A. (2023): Statistics and Geodata Analysis using R (SOGA-R). Department of Earth Sciences, Freie Universitaet Berlin. Disponível em: https://www.geo.fu-berlin.de/en/v/soga-r/Advances-statistics/Multivariate-approaches/Principal-Component-Analysis/index.html
Huang, F. L. (2019). Principal Components Analysis using R. Publicado em Oct 27, 2019. Disponível em: https://francish.net/post/pca_2020/
James, G., Witten, D., Hastie, T., Tibshirani, R. (2021). An Introduction to Statistical Learning with Applications in R Second Edition. Springer. Revisão em 2023. Disponível em: https://hastie.su.domains/ISLR2/ISLRv2_corrected_June_2023.pdf.download.html
Jolliffe, I.T., Cadima, J. (2016). Principal component analysis: a review and recent developments. Phil. Trans. R. Soc. A 374: 20150202. http://dx.doi.org/10.1098/rsta.2015.0202.
Kim, J.K. (2024). Step-by-Step Guide to Calculating and Analyzing Principal Component Analysis (PCA) by Hand. Publicado em: November 1, 2024. Disponível em: https://agronomy4future.com/archives/16454.
Margalit, D., Rabinoff, J. (2019). Interactive Linear Algebra. Publicado em June 3, 2019. Disponível em: https://textbooks.math.gatech.edu/ila/eigenvectors.html.