Análise de Componentes Principais - ACP

---

Exemplo 1 - Taxa de Delitos por 100000 habitantes por divisão territorial das polícias do Estado de SP (Deinter) em 2002.

library(factoextra)
library(FactoMineR)
library(ggplot2)
library(ggcorrplot)
library(psych)

# Deinter
deinter  <- c("SJRP", "RP", "Bauru", "Campinas", "Sorocaba", "SP", "SJC", "Santos", "GSP")
# X1
homicidio_doloso <- c(10.85,14.13,8.62,23.04,16.04,43.74,25.39,42.86,42.55)
# X2
furto <- c(1500.80,1496.07,1448.79,1277.33,1204.02,1190.94,1292.91,1590.66,797.16)
# X3
roubo <- c(149.35,187.99,130.97,424.87,214.36,1139.52,358.39,721.90,520.73)
roubo_furto_veiculos <- c(108.38,116.66,69.98,435.75,207.06,909.21,268.24,275.89,602.63)
dados <- data.frame(homicidio_doloso,furto,roubo,roubo_furto_veiculos)
row.names(dados) <- deinter
str(dados)

'data.frame':   9 obs. of  4 variables:
 $ homicidio_doloso    : num  10.85 14.13 8.62 23.04 16.04 ...
 $ furto               : num  1501 1496 1449 1277 1204 ...
 $ roubo               : num  149 188 131 425 214 ...
 $ roubo_furto_veiculos: num  108 117 70 436 207 ...

dados

#tail(dados)
# sumário estatístico 
describe(dados)

                     vars n    mean     sd  median trimmed    mad    min     max
homicidio_doloso        1 9   25.25  14.36   23.04   25.25  18.07   8.62   43.74
furto                   2 9 1310.96 239.48 1292.91 1310.96 231.11 797.16 1590.66
roubo                   3 9  427.56 330.76  358.39  427.56 252.64 130.97 1139.52
roubo_furto_veiculos    4 9  332.64 275.01  268.24  332.64 237.01  69.98  909.21
                       range  skew kurtosis     se
homicidio_doloso       35.12  0.27    -1.84   4.79
furto                 793.50 -0.82    -0.33  79.83
roubo                1008.55  0.97    -0.37 110.25
roubo_furto_veiculos  839.23  0.91    -0.56  91.67

# gráfico boxplot 
boxplot(dados) # dados originais

boxplot(scale(dados)) # dados padronizados

correlações

# Correlation matrix
matcor <- round(cor(dados), 2)
matcor

                     homicidio_doloso furto roubo roubo_furto_veiculos
homicidio_doloso                 1.00 -0.44  0.88                 0.80
furto                           -0.44  1.00 -0.26                -0.63
roubo                            0.88 -0.26  1.00                 0.88
roubo_furto_veiculos             0.80 -0.63  0.88                 1.00

# Plot
ggcorrplot(matcor, hc.order = TRUE, 
           type = "lower", 
           lab = TRUE, 
           lab_size = 3, 
           method="circle", 
           colors = c("tomato2", "white", "springgreen3"), 
           title="Correlograma", 
           ggtheme=theme_bw)

teste de Bartlett

Bartlett.sphericity.test

Bartlett.sphericity.test <- function(x)
{
  method <- "Bartlett's test of sphericity"
  data.name <- deparse(substitute(x))
  x <- subset(x, complete.cases(x)) # Omit missing values
  n <- nrow(x)
  p <- ncol(x)
  chisq <- (1-n+(2*p+5)/6)*log(det(cor(x)))
  df <- p*(p-1)/2
  p.value <- pchisq(chisq, df, lower.tail=FALSE)
  names(chisq) <- "X-squared"
  names(df) <- "df"
  return(structure(list(statistic=chisq, parameter=df, p.value=p.value,
                        method=method, data.name=data.name), class="htest"))
}
Bartlett.sphericity.test(dados)

    Bartlett's test of sphericity

data:  dados
X-squared = 33.089, df = 6, p-value = 1.008e-05

PCA

# PCA com a matriz de cov
res.pca.cov <- PCA(dados, scale.unit = F, graph = FALSE)
# matriz de covariância
round(cov(dados),2)

                     homicidio_doloso     furto     roubo roubo_furto_veiculos
homicidio_doloso               206.07  -1526.24   4190.21              3156.38
furto                        -1526.24  57352.79 -20611.99            -41428.18
roubo                         4190.21 -20611.99 109401.24             80241.60
roubo_furto_veiculos          3156.38 -41428.18  80241.60             75628.46

# autovalores
round(res.pca.cov$eig,3)

# autovetores
round(res.pca.cov$svd$V,3)

       [,1]   [,2]   [,3]   [,4]
[1,]  0.029  0.006 -0.117  0.993
[2,] -0.310  0.866  0.389  0.050
[3,]  0.716  0.484 -0.496 -0.082
[4,]  0.624 -0.125  0.768  0.073

Na primeira CP os coeficientes das variáveis roubo e roubo e furto de veículos são os maiores, pois as variâncias destas varáveis são maiores que as demais.

# A proporção de variação retida pelos componentes principais (CP) pode ser extraída da seguinte forma
round(res.pca.cov$eig,2)

# A importância dos CP pode ser visualizada usando o scree plot :
fviz_screeplot(res.pca.cov, ncp=4)+ theme_minimal()

# A correlação entre uma variável e um CP é chamada de carga (loadings). 
round(res.pca.cov$var$cor,4)

                       Dim.1   Dim.2   Dim.3   Dim.4
homicidio_doloso      0.8749  0.0956 -0.3933  0.2661
furto                -0.5624  0.8231  0.0784  0.0008
roubo                 0.9401  0.3331 -0.0723 -0.0010
roubo_furto_veiculos  0.9855 -0.1037  0.1347  0.0010

As CP também podem ser obtidas a partir as variáveis padronizadas

# PCA com a matriz de cor
res.pca.cor <- PCA(dados, scale.unit = T, graph = FALSE)
# matriz de covariância
round(cor(dados),4)

                     homicidio_doloso   furto   roubo roubo_furto_veiculos
homicidio_doloso               1.0000 -0.4440  0.8825               0.7995
furto                         -0.4440  1.0000 -0.2602              -0.6290
roubo                          0.8825 -0.2602  1.0000               0.8822
roubo_furto_veiculos           0.7995 -0.6290  0.8822               1.0000

# autovalores
round(res.pca.cor$eig,3)

# autovetores
round(res.pca.cor$svd$V,3)

       [,1]   [,2]   [,3]   [,4]
[1,]  0.533  0.213 -0.769  0.283
[2,] -0.361  0.870  0.108  0.317
[3,]  0.526  0.440  0.233 -0.690
[4,]  0.557 -0.056  0.586  0.586

# A proporção de variação retida pelos componentes principais (CP) pode ser extraída da seguinte forma
res.pca.cor$eig

# A importância dos CP pode ser visualizada usando o scree plot :
fviz_screeplot(res.pca.cor, ncp=4)+ theme_minimal()

# A correlação entre uma variável e um CP é chamada de carga (loadings). 
round(res.pca.cor$var$cor,4)

                       Dim.1   Dim.2   Dim.3   Dim.4
homicidio_doloso      0.9238  0.1904 -0.3312  0.0248
furto                -0.6252  0.7786  0.0465  0.0279
roubo                 0.9116  0.3939  0.1003 -0.0605
roubo_furto_veiculos  0.9649 -0.0501  0.2525  0.0515

O quadrado da correlação entre a variável e a CP representa a porcentagem de variância de uma das variáveis originais explicada por uma das CP

Conceito de Comunalidade

a comunalidade é a soma dos quadrados das correlações entre cada variável i e a componente principal j (ou o mesmo que o índice cos2). A soma limite-se ao número do componenentes retidos. Em nosso exemplo ilustrativo retemos todos os componentes que é igual ao número de variáveis.

# banseando-se na matriz de cov
round(res.pca.cov$var$cor^2,4)

                      Dim.1  Dim.2  Dim.3  Dim.4
homicidio_doloso     0.7654 0.0091 0.1547 0.0708
furto                0.3163 0.6775 0.0061 0.0000
roubo                0.8838 0.1110 0.0052 0.0000
roubo_furto_veiculos 0.9711 0.0107 0.0181 0.0000

# ou  round(res.pca.cov$var$cos2,4)
# veja que reter apens um CP é suficiente para explicar cerca de 76,5% da variabilidade do homicidio_doloso, e 31,6% da variável furto.
# a soma é igual a 100%
sum((res.pca.cov$var$cor^2)[1,])

[1] 1

# sum(res.pca.cov$var$cos2[1,])
# banseando-se na matriz de cor
round(res.pca.cor$var$cor^2,4)

                      Dim.1  Dim.2  Dim.3  Dim.4
homicidio_doloso     0.8534 0.0363 0.1097 0.0006
furto                0.3908 0.6062 0.0022 0.0008
roubo                0.8311 0.1552 0.0101 0.0037
roubo_furto_veiculos 0.9311 0.0025 0.0637 0.0026

# ou  round(res.pca.cor$var$cos2,4)
# veja que reter apens um CP é suficiente para explicar cerca de 85,3% da variabilidade do homicidio_doloso, e 39,1% da variável furto.

Mapa Fatorial

Quando um subespaço projetivo bidimensional determinado por duas direções principais escolhidas (CP), sua imagem geométrica plana com os pontos projetados e o círculo de correlações é denominada MAPA FATORIAL

# Quanto mais próxima uma variável for do círculo de correlações, melhor sua representação no mapa fatorial (e mais importante é a variável para a interpretação desses componentes)
# As variáveis próximas ao centro do gráfico são menos importantes para os primeiros componentes.
# No gráfico abaixo os componentes são coloridas de acordo com os valores do coseno quadrado:
fviz_pca_var(res.pca.cor, col.var="cos2") +
scale_color_gradient2(low="white", mid="blue", 
                    high="red", midpoint=0.5) + theme_minimal()

# Coordenadas de variáveis
round(res.pca.cor$var$coord,2)

                     Dim.1 Dim.2 Dim.3 Dim.4
homicidio_doloso      0.92  0.19 -0.33  0.02
furto                -0.63  0.78  0.05  0.03
roubo                 0.91  0.39  0.10 -0.06
roubo_furto_veiculos  0.96 -0.05  0.25  0.05

# Cos2: é uma medida que indica a qualidade da representação para variáveis no mapa fatorial
round(res.pca.cor$var$cos2,2)

                     Dim.1 Dim.2 Dim.3 Dim.4
homicidio_doloso      0.85  0.04  0.11     0
furto                 0.39  0.61  0.00     0
roubo                 0.83  0.16  0.01     0
roubo_furto_veiculos  0.93  0.00  0.06     0

Contribuições das variáveis para os componentes principais

As variáveis que são correlacionadas com PC1 e PC2 são as mais importantes para explicar a variabilidade no conjunto de dados. Variáveis que não se correlacionam com nenhum PC ou correlacionadas com as últimas dimensões são variáveis com baixa contribuição e podem ser removidas para simplificar a análise geral. As contribuições das variáveis na contabilização da variabilidade em uma determinada componente principal são (em porcentagem): (variável.cos2 * 100) / (cos2 total da componente)

# A contribuição das variáveis pode ser extraída da seguinte forma:
round(res.pca.cor$var$contrib,2)

                     Dim.1 Dim.2 Dim.3 Dim.4
homicidio_doloso     28.39  4.53 59.08  8.00
furto                13.00 75.76  1.17 10.07
roubo                27.64 19.39  5.42 47.55
roubo_furto_veiculos 30.97  0.31 34.34 34.38

# veja que a soma é igual a 100%
sum (res.pca.cor$var$contrib[,1])  

[1] 100

# Quanto maior o valor da contribuição, mais a variável contribui para o componente.
# As variáveis mais importantes associadas a um determinado PC podem ser visualizadas, usando a função fviz_contrib () [factoextra package], da seguinte forma:
# Contribuições de variáveis no PC1
fviz_contrib(res.pca.cor, choice = "var", axes = 1)+ theme_minimal()

# Contribuições de variáveis no PC2
fviz_contrib(res.pca.cor, choice = "var", axes = 2)+ theme_minimal()

# Contribuição total nos PC1 e PC2
fviz_contrib(res.pca.cor, choice = "var", axes = 1:2)+ theme_minimal()

# Controle as cores das variáveis usando suas contribuições
# a cor representa a contribuição conjunta dim1-dim2
fviz_pca_var(res.pca.cor, col.var="contrib")+ theme_minimal()

# Alterar a cor 
fviz_pca_var(res.pca.cor, col.var="contrib") +
scale_color_gradient2(low="white", mid="blue", 
                  high="red", midpoint=50) + theme_minimal()

# A função dimdesc () [em FactoMineR] pode ser usada para identificar as variáveis mais correlacionadas com uma determinada componente principal.
res.desc <- dimdesc(res.pca.cor, axes = c(1,2))
# Descrição da dimensão 1
res.desc$Dim.1

$quanti
                     correlation      p.value
roubo_furto_veiculos   0.9649332 2.569051e-05
homicidio_doloso       0.9238085 3.729002e-04
roubo                  0.9116461 6.186133e-04

# Descrição da dimensão 2
res.desc$Dim.2

$quanti
      correlation    p.value
furto   0.7786169 0.01343302

# Descrição da dimensão 3
res.desc$Dim.3

NULL

Se a análise for baseada na matriz de correlação, a primeira componente pode ser interpretada como a média de homicídio doloso, roubo e roubo e furto de veículos, os três crimes que incluem a presença da vítima na ocorrência do mesmo. A Segunda componente principal seria interpretada como a variável furto.

Análise de pontos (escores, objetos, indivíduos)

# Gráfico de escores (indivíduos ou pontos objetos)
# As coordenadas dos escores nos componentes principais são:
round(res.pca.cor$ind$coord,2)

         Dim.1 Dim.2 Dim.3 Dim.4
SJRP     -1.82  0.16  0.19  0.07
RP       -1.60  0.25  0.05  0.07
Bauru    -1.94 -0.09  0.20 -0.09
Campinas  0.18 -0.19  0.34  0.15
Sorocaba -0.82 -0.83  0.03 -0.16
SP        3.36  0.71  0.73 -0.05
SJC      -0.22 -0.15 -0.21 -0.02
Santos    0.62  1.78 -0.78 -0.02
GSP       2.24 -1.64 -0.55  0.04

fviz_pca_ind(res.pca.cor)+ theme_minimal()

#Cos2: qualidade da representação para escores nos componentes principais
# O coseno quadrado mostra a importância de um componente para uma determinada observação.
round(res.pca.cor$ind$cos2,3)

         Dim.1 Dim.2 Dim.3 Dim.4
SJRP     0.980 0.008 0.011 0.002
RP       0.974 0.023 0.001 0.002
Bauru    0.986 0.002 0.010 0.002
Campinas 0.163 0.175 0.559 0.103
Sorocaba 0.485 0.497 0.001 0.017
SP       0.916 0.041 0.043 0.000
SJC      0.414 0.195 0.390 0.002
Santos   0.093 0.763 0.144 0.000
GSP      0.627 0.335 0.038 0.000

fviz_pca_ind(res.pca.cor, col.ind="cos2") +
scale_color_gradient2(low="white", mid="blue", 
   high="red", midpoint=0.50) + theme_minimal()

# Contribuição dos escores para os componentes principais 
round(res.pca.cor$ind$contrib,2)

         Dim.1 Dim.2 Dim.3 Dim.4
SJRP     12.25  0.36  2.25  7.91
RP        9.47  0.85  0.17  6.91
Bauru    13.89  0.12  2.28 12.17
Campinas  0.12  0.50  6.93 30.64
Sorocaba  2.49  9.59  0.05 34.79
SP       41.66  6.98 31.60  4.01
SJC       0.18  0.32  2.74  0.33
Santos    1.43 44.15 35.97  0.46
GSP      18.51 37.13 18.02  2.77

# Contribuições de escores para PC1
fviz_contrib(res.pca.cor, choice = "ind", axes = 1)+ theme_minimal()

# Contribuições de escores para PC2
fviz_contrib(res.pca.cor, choice = "ind", axes = 2)+ theme_minimal()

# Contribuição total em PC1 e PC2
fviz_contrib(res.pca.cor, choice = "ind", axes = 1:2)+ theme_minimal()

# Contribuições dos escores para PC1  (apenas os "top")
fviz_contrib(res.pca.cor, choice = "ind", axes = 1:2, top = 5)+ theme_minimal()

# Mundando a cor
fviz_pca_ind(res.pca.cor, col.ind="contrib") +
scale_color_gradient2(low="white", mid="blue", 
                  high="red", midpoint=50) + theme_minimal()

O gráfico de dispersão dos escores dos dois primeiros componentes baseados na matriz de correlações juntamente com os respectivos autovetores

Este gráfico é chamado de biplot. É uma representação bidimensional de dados multivariados.

fviz_pca_biplot(res.pca.cor) + theme_minimal()

sumário

# sumário
facto_summarize(res.pca.cor, "var") # para variáveis

facto_summarize(res.pca.cor, "ind") # para escores