Introdução

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.

Diretório de trabalho

Definimos o diretório de trabalho:

setwd("/home/gf/Scripts/Moyses")

Bibliotecas

Carregamos as bibliotecas necessárias:

Importação dos dados

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

Inspeção de dados

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

Teste de Bartlett

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.

Padronização dos dados

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)

Passo 1 - Obtenção da matriz de covariância amostral

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

Passo 2. Obtenção dos autovetores e autovalores

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

Passo 3. Percentual de explicação acumulada de cada componente principal

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
  }
}

Visualização dos resultados

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.

Diagrama de dispersão

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)

Verificação da ortogonalidade dos componentes

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

correlação entre as variáveis originais e os componentes principais

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.