Análise de Componentes Principais (1)
Prof. Lupercio F. Bessegato
10/10/2020
Análise de Componentes Principais - Cálculo Matricial
Inicialmente serão apresentados os cálculos matriciais para todas as quantidades e gráficos envolvidos com a análise de componentes principais. Recomend-ase fortemente que todos entendam as operações matriciais apresentadas e que consigam utilizá-las para realizar a análise de componentes prinicpais sem necessidade de uso de comandos ou pacotes específicos. A compreensão da metodologia auxiliará no entendimento de análise fatorial, que envolve modelgame probabilística.
Conjunto de dados iris:
Conjunto de dados buil-in do R base, com 150 observações, com 4 variáveis relacionadas com dimensões morfólogicas de flores íris, com 50 observações de três de suas espécies: "setosa", "versicolor" e `“virginica”’. O vetor dde médias (\(\bar{\mathbf{x}}\)) e a matriz de covariâncias (\(\mathbf{S}\)) amostrais das quatro variáveis do conjunto de dados são:
# vetor de médias
x.bar <- sapply(iris[-5], mean)
x.bar## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 5.843333 3.057333 3.758000 1.199333 # matriz de covariâncias
S <- cov(iris[-5])
S## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 0.6856935 -0.0424340 1.2743154 0.5162707
## Sepal.Width -0.0424340 0.1899794 -0.3296564 -0.1216394
## Petal.Length 1.2743154 -0.3296564 3.1162779 1.2956094
## Petal.Width 0.5162707 -0.1216394 1.2956094 0.5810063Percebe-se que a variância da variável Petal.Length é bem maior em relação às demais e ela pode dominar a análise de componentes principais pela matriz de covariâncias dos dados.
As componentes principais das quatro variáveis métricas podem ser calculadas com o comando eigen() que fornece os autovalores e autovetores da matriz de covariâncias ou de correlações. Se consideramos a matriz de covariâncias \(\mathbf{S}\), a estrutura de autovalores e autovetores pode ser calculada, sendo utilizada para estimar as componentes principais.
iris.av <- eigen(S)
iris.av$values## [1] 4.22824171 0.24267075 0.07820950 0.02383509iris.av$vectors## [,1] [,2] [,3] [,4]
## [1,] 0.36138659 -0.65658877 -0.58202985 0.3154872
## [2,] -0.08452251 -0.73016143 0.59791083 -0.3197231
## [3,] 0.85667061 0.17337266 0.07623608 -0.4798390
## [4,] 0.35828920 0.07548102 0.54583143 0.7536574A média de cada componente é facilmente calculada pela expressão \(\hat{\mathbf{E}}^\prime \bar{\mathbf{x}}\), em que \(\hat{\mathbf{E}}_{4\mathrm{x}4}\) é a matriz dos autovetores amostrais, ou seja, matriz cujas coluna são os quatro autovetores da decomposição espectral da matriz \(\mathbf{S}\) e \(\bar{\mathbf{x}}\) é o vetor de médias amostrais do conjunto de dados em análise. Assim, o vetor de médis amostral das componentes amostrais é:
# média das componentes principais amostrais
escores.media <- t(iris.av$vectors) %*% x.bar
escores.media## [,1]
## [1,] 5.50236513
## [2,] -5.32695258
## [3,] -0.63185272
## [4,] -0.03335171Por sua vez, a variância total da matriz \(S\) é 4.573 sendo distribuída por cada componente de acordo com seus autovalores. Assim, a porcentagem da variância totl explicada por cada componente é:
iris.av$values/sum(iris.av$values)## [1] 0.924618723 0.053066483 0.017102610 0.005212184indicando que que \(92\%\) da variabilidade total é explicada pela primeira componente. O scree plot do ajuste é obtido pelo comando abaixo:
# Scree plot
plot(iris.av$values, type = "b", pch = 21, bg = "black")
Por outro lado, da observação dos elementos do autovetores (coeficientes das componentes principais) percebe-se entretanto que a variável Sepal.Widthcarrega pouco na primeira componente, estando mais presente na segunda. Um gráficos dos coeficientes das componentes auxiliam na visualização desses coeficientes e na interpretação das componentes:
plot(iris.av$vectors[, 1:2], pch = 21, bg = "black", xlim = c(-0.2, 1.0),
ylim = c(-0.8, 0.3), xlab = "Componente 1", ylab = "Componente 2",
main = "Coeficientes das Variáveis")
abline(h = 0, v = 0, lty = 1)
arrows(0, 0, iris.av$vectors[, 1], iris.av$vectors[, 2], length = 0.25,
angle = 15, lty = 3)
text(iris.av$vectors[,1:2], labels = names(iris[-5]), pos = c(1, 1, 3, 3),
cex = 0.85)
A correlação entre a i-ésima componente e k-ésima variável é obtida pela expressão abaixo: \[ \begin{equation} \tag{1} r(Y_i, X_k) = \frac{\hat{e}_{ik}\sqrt{\hat{\lambda}_i}}{\sqrt{s_{kk}}} . \end{equation} \]
Consideremos as seguintes matrizes: \(\hat{\mathbf{E}}\), a matriz dos autovetores amostrais, ou seja, aquela cujas colunas são os autovetores de \(\mathbf{S}\), \(\hat{\boldsymbol{\Lambda}}\), a matriz diagonal dos autovalores amostrais e \(\mathbf{V}\), a matriz diagonal das variâncias amostrais de cada variável (diagonal de \(\mathbf{S}\)). Assim, a expressão \(\tag{1}\) pode ser obtida matricialmente, forecendo todas as \(4^2\) correlações amostrais entre as componentes e as variáveis originais; \[ \mathbf{V}^{-1/2}\, \hat{\mathbf{E}} \, \hat{\boldsymbol{\Lambda}}^{1/2} \] Verifique essa expressão matricial e acostume-se a entender, a usar e a elaborar esse tipo de expressão (1) podem ser obtidas pelos comandos abaixo:
# matriz dos autovetores
E <- iris.av$vectors
colnames(E) <- paste("Comp.", 1:4)
rownames(E) <- names(iris[-5])
# matriz diagonal dos autovalores
Lambda <- diag(iris.av$values)
# matriz diagonal dos desvios-padrão
V.meio <- diag(1/sqrt(diag(S)))
# matriz das correlações entre as componentes e as variáveis
cor.YX <- V.meio %*% E %*% sqrt(Lambda)
colnames(cor.YX) <- paste("Comp.", 1:4)
rownames(cor.YX) <- names(iris[-5])
round(cor.YX, 3)## Comp. 1 Comp. 2 Comp. 3 Comp. 4
## Sepal.Length 0.897 -0.391 -0.197 0.059
## Sepal.Width -0.399 -0.825 0.384 -0.113
## Petal.Length 0.998 0.048 0.012 -0.042
## Petal.Width 0.967 0.049 0.200 0.153Os valores confirmam a análise dos coeficientes dos autovetores e mostram a alta correlação entre a 1ª componente e a variável Petal.Length, que domina a variância total (68.1% da variância total).
Dessa maneira, ao invés das variáveis originais, podemos trabalhar com os escores das componentes, ou sejam o valor de cada componente para todas as observações. Denotando a matriz de dados por \(\mathbf{X}_{150 \times 4}\), a matriz dos escores amostrais, \(\hat{\mathbf{Y}}_{150 \times 4}\) é obtida pela expressão matricial \(\mathbf{X} \, \hat{\mathbf{E}}\), em que \(\hat{\mathbf{E}}\) é a matriz dos autovetores amostrais de \(\mathbf{S}\). Os \(150\) escores de cada uma das quatro componentes amostrais são obtidos pelos comandos abaixo:
# Escores das observações
escores <- as.matrix(iris[-5])%*% iris.av$vectors
dim(escores)## [1] 150 4head(escores)## [,1] [,2] [,3] [,4]
## [1,] 2.818240 -5.646350 -0.6597675 -0.03108928
## [2,] 2.788223 -5.149951 -0.8423170 0.06567484
## [3,] 2.613375 -5.182003 -0.6139525 -0.01338332
## [4,] 2.757022 -5.008654 -0.6002933 -0.10892753
## [5,] 2.773649 -5.653707 -0.5417735 -0.09461031
## [6,] 3.221505 -6.068283 -0.4631751 -0.05755257As médias dos escores das componentes amostrais são as quatro combinações lineares das médias amostrais de cada variável, cujos coeficientes são os elementos dos quatro autovetores. Por outro lado, os devios padrão amostrais dos escores são a raiz quadrada de cada um dos quatro autovalores amostrais. Apenas para conferir, acompanhe os cálculos abaixo:
# medias e desvios-padrão dos escores
apply(escores, 2, function(x) c(media = mean(x), d.pad = sd(x)))## [,1] [,2] [,3] [,4]
## media 5.502365 -5.3269526 -0.6318527 -0.03335171
## d.pad 2.056269 0.4926162 0.2796596 0.15438618 # vetor de médias das componentes.
escores.media <- as.vector(t(iris.av$vectors) %*% x.bar)
escores.media## [1] 5.50236513 -5.32695258 -0.63185272 -0.03335171 # desvio-padrão das componentes
sqrt(iris.av$values)## [1] 2.0562689 0.4926162 0.2796596 0.1543862Em geral, usamos O diagrama de dispersão dos escores das duas primeiras componentes para tentarmos obter uma visualização de todas as variáveis. Como sabemos que os dados estão estratificados pelo fator Species, identificaremos cada espécie de acordo com seus níveis:
# configuração cores
library(RColorBrewer)
# display.brewer.all()
cores <- brewer.pal(n = 3, name = "Set2")
cor_species <- cores[1:3]
# plot
plot(escores[ , 1:2], col = cores[as.numeric(iris[, 5])], pch = 16,
xlab = "Componente 1", ylab = "Componente 2")
abline(v = escores.media[1], h = escores.media[2], lty = 2, col = "gray")
legend("topright", levels(iris[, 5]), pch = 16, col = cores, cex = 0.90,
bty = "n")
Ao plotarmos as observações neste sistema de eixo coordenados, agora considerando índices das dimensões morfológicas observadas dessas plantas, percebemos com mais clareza os clusters relacionados a cada espécie, já que conseguimos praticamente esboçar a região de cada uma das três espécies, neste plano cartesiano das duas primeiras componentes.
Há comandos que fornecem os escores amostrais centralizados, ou seja, todos com média zero. Os gráficos obtidos com esses escores são tem o mesmo formato, estando centralizados nas médias amostrais das componentes. Centralizamos os escores com o comando scale() e, caso tenhamos os escores centralizados, conseguimos recuperá-los com o auxílio do comando sweep.
# plot das duas primeiras componentes centralizadas
# centralização os escores
escores.centro <- apply(escores, 2,
function(x) scale(x, center = TRUE, scale = FALSE))
# plot dos escores centralizados
plot(escores.centro[ , 1:2], col = cores[as.numeric(iris[, 5])], pch = 16,
xlab = "Componente 1", ylab = "Componente 2")
abline(v = 0, h = 0, lty = 2, col = "gray")
legend("topright", levels(iris[, 5]), pch = 16, col = cores, cex = 0.90,
bty = "n")
O gráfico é o mesmo, com exceção de sua traslação para o centróide das duas primeiras componentes amostrais.
É importante salientar que o gráfico das variáveis originais Petal.Length e Petal.Width discriminam razoavelmente essas três espécies. Dêem uma olhada em seu gráfico. Qual discrimina melhor? Lembrando ainda que a variável Petal.Width contém uma grande parcela da variabilidade total.
# plot das duas variáveis originais com maior variabilidade
plot(Petal.Length ~ Petal.Width, data = iris, pch = 16,
col = cores[as.numeric(iris[, 5])])
legend("topleft", levels(iris[, 5]), pch = 16, col = cores, cex = 0.90,
bty = "n")
Desse gráfico, percebemos que, quanto à discriminação das espécies, a própria variável Petal.Widthjá auxiia bastante nesta tarefa.
Podemos então tentar visualizar as variáveis usando apenas a primeira compenente principal, que contèm \(92\%\) da variabilidade total. Aqui também estratificamos os escores da primeira componente de acordo com os níveis do fator Species.
# lista com vetores dos escores e Petal.Width por nível de Species
estrato.cp <- split(x = escores[, 1], f = iris$Species)
estrato.PW <- split(x = iris$Petal.Width, f = iris$Species)
par(mfrow = c(2,1), mar = c(2.1, 4.1, 4.1, 2.1))
# gráfico da primeira componente
stripchart(estrato.cp, main = "Primeira Componente Principal - iris",
xlab = "", ylab = "Espécie",
col = cores, pch = 16, method = "jitter",
cex.axis = 0.75)
# Gráfico de Petal.Width
stripchart(estrato.PW, main = "Petal.Width - iris",
xlab = "", ylab = "Espécie",
col = cores, pch = 16, method = "jitter",
cex.axis = 0.75)
Aparentemente, em termos de discriminação, a primeira componente discrimina de maneira similar à variável Petal.Width, que tem maior variância dentre as quatro variáveis originais do conjunto de dados.
Entretanto, para analisar a variabilidade do conjunto de dados. representada pela matriz de covariâncias amostral, podemos substituir as variáveis originais pelas componentes principais e obteremos a mesma matriz \(\mathbf{S}\), ou seja, \(\mathbf{S} = \hat{\mathbf{E}} \, \hat{\boldsymbol{\Lambda}} \, \hat{\mathbf{E}} ^\prime\). Assim:
# Reconsituição de S
# matriz dos autovetores
E <- iris.av$vectors
colnames(E) <- paste("Comp.", 1:4)
rownames(E) <- names(iris[-5])
# matriz dos autovalores
Lambda <- diag(iris.av$values)
# reconstituição de S
S.rec <- E %*% Lambda %*% t(E)
# Diferença entre as matrizes
S.rec - S## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length -3.330669e-16 9.020562e-17 -6.661338e-16 -2.220446e-16
## Sepal.Width 1.040834e-16 -8.049117e-16 1.776357e-15 1.054712e-15
## Petal.Length -6.661338e-16 1.720846e-15 -3.108624e-15 -1.554312e-15
## Petal.Width -2.220446e-16 1.054712e-15 -1.554312e-15 -7.771561e-16 # soma dos quadrados das diferenças
sum((S.rec - S)^2)## [1] 2.520474e-29As matrizes são iguais, como era de se esperar! Importante salientar que essa soma de quadrados cobre todos os elementos da matriz e não só os elementos da diagonal prinicpal, como é o caso quando focamos a variância total. Se trabalhamos apenas com uma \(k\) componentse principis, \(k <p\) teremos uma aproximação de \(\mathbf{S}\), ou seja \[\mathbf{S}_{p\times p} \approx \sum_{i=1}^k \hat{\lambda}_{i} \hat{\mathbf{e}}_i \hat{\mathbf{e}}_i^\prime = \hat{\mathbf{E}}^*_{k \times k} \hat{\boldsymbol{\Lambda}}^*_{k \times k} {\hat{\mathbf{E}}^*}^\prime_{k \times k}. \]
No exemplo, a primeira componente principal contém \(92\%\) da variância total dos dados, assim, ao utilizarmos apenas ela teremos:
# matrizes reduzidas da decomposição
E.k <- as.matrix(E[, 1])
Lambda.k <- Lambda[1, 1]
# aproximação de S
S.k <- E.k %*% Lambda.k %*% t(E.k)
# diferença entre matrizes
S.k - S## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length -0.13348401 -0.08671892 0.034702829 0.031205898
## Sepal.Width -0.08671892 -0.15977263 0.023498048 -0.006406600
## Petal.Length 0.03470283 0.02349805 -0.013236687 0.002189455
## Petal.Width 0.03120590 -0.00640660 0.002189455 -0.038222019 # soma dos quadrados das diferenças
sum((S.k - S)^2)## [1] 0.06557393Uma maneira mais simples de entendermo esse processo de decomposição/reconstituição da matriz \(\mathbf{S}\) é estruturarmos a reconstituição de \(\mathbf{S}\), com a expressão abaixo. A matriz \(\hat{\boldsymbol{\Lambda}}^{1/2}\) é a matriz diagonal das raízes quadrados dos autovalores de \(\mathbf{S}\): \[\begin{align} \mathbf{S} = \hat{\mathbf{E}} \, \hat{\boldsymbol{\Lambda}} \, \hat{\mathbf{E}} ^\prime = \left(\hat{\mathbf{E}} \, \hat{\boldsymbol{\Lambda}}^{1/2} \right) \left(\hat{\boldsymbol{\Lambda}}^{1/2} \hat{\mathbf{E}}^\prime \right) = \left(\hat{\mathbf{E}} \, \hat{\boldsymbol{\Lambda}}^{1/2} \right) \left(\hat{\mathbf{E}} \, \hat{\boldsymbol{\Lambda}}^{1/2} \right)^\prime = \mathbf{L} \, \mathbf{L}^\prime. \end{align}\]
A matriz \(\mathbf{L} = \hat{\mathbf{E}} \, \hat{\boldsymbol{\Lambda}}^{1/2}\) é denominada por alguns autores matriz de loadings. Notem que seus elementos são os numeradores das correlações entre as variáveis originais e as componentes, ou seja, \(\hat{e}_{ik} \sqrt{\hat{\lambda}_i}\), para a i-ésima componente prinicpal e k-ésima variável. Ela será importante nas verificações dos ajustes e, em algumas ocasiões tem de ser calculada, pois, alguns comandos não a fornece.
Estruturando com a matrz de loadings, Os cálculos ficam mais simplificados. Vejam a repetição do procedimento anterior para a matriz de covariâncias \(\mathbf{S}\):
# -------------------------------
# Matriz de loadings
Lambda.meio <- diag(sqrt(iris.av$values))
# matriz de loadings
L <- E %*% Lambda.meio
colnames(L) <- paste("Comp.", 1:4)
rownames(L) <- names(iris[-5])
L## Comp. 1 Comp. 2 Comp. 3 Comp. 4
## Sepal.Length 0.7431080 -0.32344628 -0.16277024 0.04870686
## Sepal.Width -0.1738010 -0.35968937 0.16721151 -0.04936083
## Petal.Length 1.7615451 0.08540619 0.02132015 -0.07408051
## Petal.Width 0.7367389 0.03718318 0.15264701 0.11635429 # reconstrução de S
(S.rec <- L %*% t(L))## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 0.6856935 -0.0424340 1.2743154 0.5162707
## Sepal.Width -0.0424340 0.1899794 -0.3296564 -0.1216394
## Petal.Length 1.2743154 -0.3296564 3.1162779 1.2956094
## Petal.Width 0.5162707 -0.1216394 1.2956094 0.5810063 # variância total
sum(diag(S.rec))## [1] 4.572957sum(diag(S))## [1] 4.572957 # soma dos quadrados das diferenças entre as matrizes
sum((S.rec - S)^2)## [1] 3.48941e-29Da mesma maneira, pode-se aproximar a matriz de covariâncias amostral \(\mathbf{S}\) para uma quantidade \(k\), \(k < p\) componentes principais, usando uma matriz composta pelas \(k\) primeiras coluna da matriz de loadings \(\mathbf{L}\), de maneira que \[ \mathbf{S} \approx \mathbf{L}^*_{p\times k} \, {\mathbf{L}^*}^\prime_{k \times p} \, , \]
em que \(\mathbf{L}^*_{p\times k}\) é a matriz composta pelas k primeiras colunas de \(\mathbf{L}\). Dessa maneira, a aproximação de \(\mathbf{S}\) com a primeira componente principal é calculada a seguir, alcançando os mesmos resultados obtidos anteriormente:
# Aproximação de S com CP1
(S.1 <- L[, 1] %*% t(L[, 1]))## Sepal.Length Sepal.Width Petal.Length Petal.Width
## [1,] 0.5522095 -0.12915293 1.3090183 0.5474766
## [2,] -0.1291529 0.03020679 -0.3061583 -0.1280460
## [3,] 1.3090183 -0.30615833 3.1030412 1.2977989
## [4,] 0.5474766 -0.12804597 1.2977989 0.5427842 # variância total aproximada
sum(diag(S.1))## [1] 4.228242sum(diag(S))## [1] 4.572957 # % variância total explicada por CP1
sum(diag(S.1))/sum(diag(S))## [1] 0.9246187 # Soma dos quadrados da diferenças dos elementos de S
sum((S.1 - S)^2)## [1] 0.06557393Pode-se verificar assim a contribuição das \(k\) primeiras componentes na reconstituição da variância de cada uma das variáveis. A diagonal da matriz \(\mathbf{S}^*\) informa a quantidade da variância de cada variável que foi recuperada pelas \(k\) componentes. Essas quantidades são denominadas comunalidades. No nosso exemplo, \(k = 1\) e as comunalidades são:
# comunalidades
comunalidade <- diag(S.1)
names(comunalidade) <- names(iris[-5])
comunalidade## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.55220950 0.03020679 3.10304116 0.54278425E a porcentagem de variância de cada variável explicada pela Componente 1 é:
# % variância de cada variável explicada por CP1
explicada <- comunalidade/diag(S)
explicada## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.8053299 0.1590003 0.9957524 0.9342141ou seja, a componente capturou praticamente toda a variância de Petal.Length, uma quantidade muito grande das variáveis Petal.Width e Sepal.Length, mas apenas 15.9% da variância da variável Sepal.Width. Poderíamos assim calcular uma variância “específica” em cada variável, de maneira que sua variânica pudesse ser dividida em uma quantidade que é reconstruída (explicada) pelas componentes principais e uma paracela que não é comum às componentes, sendo considerada “específica” da variável (não caprutada pelas componentes). Importante ressaltar que se forem usadas todas as \(p\) componentes as variâncias específicas serão nulas. Em nosso exemplo, com apenas a primeira componente principal, as variâncias específicas de cada variável serão:
# variância "específica"
diag(S - S.1)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.13348401 0.15977263 0.01323669 0.03822202 # proproção específica de cada variável
diag(S - S.1)/diag(S)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.194670078 0.840999656 0.004247595 0.065785898Repetiremos o procedimento completo para o caso de usarmos as duas primeiras componentes principais:
# ----------------------------------------
# Aproximação de S com CP1 e CP2
(S.2 <- L[, 1:2] %*% t(L[, 1:2]))## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 0.65682700 -0.01281273 1.281394 0.5354498
## Sepal.Width -0.01281273 0.15958324 -0.336878 -0.1414204
## Petal.Length 1.28139395 -0.33687803 3.110335 1.3009745
## Petal.Width 0.53544983 -0.14142037 1.300975 0.5441668 # variância total aproximada
sum(diag(S.2))## [1] 4.470912sum(diag(S))## [1] 4.572957 # % variância total explicada por CP1
sum(diag(S.2))/sum(diag(S))## [1] 0.9776852 # Soma dos quadrados da diferenças dos elementos de S
sum((S.2 - S)^2)## [1] 0.006684838 # comunalidades
comunalidade <- diag(S.2)
names(comunalidade) <- names(iris[-5])
comunalidade## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.6568270 0.1595832 3.1103354 0.5441668 # % variância de cada variável explicada por CP1
explicada <- comunalidade/diag(S)
explicada## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.9579017 0.8400028 0.9980931 0.9365937 # variância "específica"
diag(S - S.2)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.028866511 0.030396181 0.005942471 0.036839430 # proproção específica de cada variável
diag(S - S.2)/diag(S)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.042098270 0.159997233 0.001906913 0.063406253A soma dos quadrados das diferenças dos elementos da matriz de covariâncias \(\mathbf{S}\) e da matriz aproximada pelas duas primeiras componentes principais diminuiu para 0.0066848, quando perfazia 0.0655739 para a matriz aproximada apemnas pela primeira componente. Além disso, a porcentagem das comunalidades passou de 80% em todas as variáveis, diminuindo as variâncias específicas, ou sejam a quantidade de variabilidade que não foram explicadas pelas duas componentes. De maneira geral, a introdução da segunda componente melhorou a proximação de \(\mathbf{S}\), devido principalmente à incorporação de uma quantidade bem maior da variabilidade da variável Sepal.Width. Note também que a matriz \(\mathbf{L}^*\) traz em seus elementos as informações sobre as comunalidades e sobre como cada uma das duas componentes colaborou na reconstituição da variância de cada uma das variáveis originais. Por exemplo, a comunalidade da variável Sepal.Width pode ser obtida com a soma dos quadrados dos elementos da segunda linha da matriz de loadings, ou seja:
sum(L[2, 1:2]^2)## [1] 0.1595832Poderíamos assim, conhecer qual a quantidade de variância cada uma das duas componentes principais contribuiu na reconstiruição da variância de cada variável original. A reconstituição será sempre parcial, a menos que sejam utilizadas todas as componentes principais de \(\mathbf{S}\). Observe abaixo o detalhamento da reconstrução de \(\mathbf{S}\), quando utilizamos as duas primeiras componentes principais:
# componentes das comunalidades
Soma <- apply(L[, 1:2]^2, 1, sum)
Soma ## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.6568270 0.1595832 3.1103354 0.5441668 # quadro da reconstrução de S
cbind(L[, 1:2]^2, Comunalidades = Soma, Especifica = diag(S) - Soma, Variancia = diag(S))## Comp. 1 Comp. 2 Comunalidades Especifica Variancia
## Sepal.Length 0.55220950 0.104617498 0.6568270 0.028866511 0.6856935
## Sepal.Width 0.03020679 0.129376444 0.1595832 0.030396181 0.1899794
## Petal.Length 3.10304116 0.007294217 3.1103354 0.005942471 3.1162779
## Petal.Width 0.54278425 0.001382589 0.5441668 0.036839430 0.5810063Nesse processo de reconstituição da matriz de covariâncias, percebe-se que 10 % da variabilidade total recuperada veio da carga da variável Petal.Length na Componente 1.
A matriz de loadings, \(\mathbf{L}\) facilita também o cálculo das correlações amostrais entre as componentes principais e as variáveis. O cálculo da expressão \((1)\) para cada uma das 16 correlações é calculado matricialmente como se segue, apresentando o mesmo resultado apresentado anteriormente:
# -------------------------------------
# correlação entre as variáveis e as componentes
# Matriz diagonal do inverso dos desvios padrão das variáveis
V.meio <- diag(1/sqrt(diag(S)))
# correlações entre Yi e Xk
cor.YX <- V.meio %*% L
rownames(cor.YX) <- rownames(S)
cor.YX## Comp. 1 Comp. 2 Comp. 3 Comp. 4
## Sepal.Length 0.8974018 -0.3906044 -0.19656672 0.05882002
## Sepal.Width -0.3987485 -0.8252287 0.38363030 -0.11324764
## Petal.Length 0.9978739 0.0483806 0.01207737 -0.04196487
## Petal.Width 0.9665475 0.0487816 0.20026170 0.15264831round(cor.YX, 3)## Comp. 1 Comp. 2 Comp. 3 Comp. 4
## Sepal.Length 0.897 -0.391 -0.197 0.059
## Sepal.Width -0.399 -0.825 0.384 -0.113
## Petal.Length 0.998 0.048 0.012 -0.042
## Petal.Width 0.967 0.049 0.200 0.153Uma outra verificação facilita a análise dos resultados. As colunas da matriz dos coeficientes das componentes são os autovetores de \(\mathbf{S}\). Assim, essa matriz é ortogonal (é um matriz especial, quais são suas propriedades?) e cada coluna tem comprimento unitário. Utilize o código abaixo e verifique os porcentuais de contribuição de cada variável nas componentes principais. O comando apply(contribuicao, 2, sum) soma os quadrados dos coeficientes de cada elemento dos autovetores, esse total será um em todas as colunas. Verifique.
contribuicao <- E**2
tabela <- rbind(contribuicao, Total = apply(contribuicao, 2, sum))
tabelaVocê perceberá que a variável Petal.Lengthcontribuiu muito com a primeira componente (73.4%) e, por sua vez, a variável Sepal.Width, muito pouco (0.7%). Podemos construir gráficos de barras dessas contribuições, ordenadas da maior para a menor, como por exemplo os gráficos a seguir, relacionados com as contribuições das variáveis para cada uma das duas primeiras componentes:
# Gráfico das Contribuições às Componentes
contribuicao <- E**2
# ---------- gráficos contribuicoes
par(mfrow=c(1,2), oma = c(0, 0, 2, 0))
# CP1
# ordenação das barras por contribuição
ordem.CP1 <- order(contribuicao[, 1], decreasing = TRUE)
# gráfico de barras das contribuições de CP1
bp1 <- barplot(contribuicao[ordem.CP1, 1]*100, beside = T, xaxt = "n",
ylab = "Contribuições (%)", ylim = c(0, 75),
main = "", sub = "CP1")
# nome das variáveis ordenadas
labels <- rownames(contribuicao)[ordem.CP1]
# abreviação das variáveis
caps <- gsub("[^A-Z]","", labels)
# nome das variáveis
text(cex = 0.85, x = bp1, y = -1.5, caps, xpd = TRUE, srt = 45, pos = 2)
# marca dos 10% de contribuição
abline(h = 10, lty = 2)
# CP2
# ordenação das barras por contribuição
ordem.CP2 <- order(contribuicao[, 2], decreasing = TRUE)
# gráfico de barras das contribuições de CP1
bp2 <- barplot(contribuicao[ordem.CP2, 2]*100, beside = T, xaxt = "n",
ylab = "Contribuições (%)", ylim = c(0, 75),
main = "", sub = "CP2")
# nome das variáveis ordenadas
labels <- rownames(contribuicao)[ordem.CP2]
# abrevisação das variáveis
caps <- gsub("[^A-Z]","", labels)
# nome das variáveis
text(cex = 0.85, x = bp2, y = -1.5, caps, xpd = TRUE, srt = 45, pos = 2)
# marca dos 10% de contribuição
abline(h = 10, lty = 2)
# Título do gráfico
mtext("Contribuições das Variáveis às Componentes", side = 3, line = 0,
outer = T, font = 2, cex = 1.5)
Os 10% é meramente indicativo. Alguns autores recomendam se descartar da interpretação das componentes as variáveis com contribuição menor que esse limiar. Percebe-se no gráfico a predominância da variável Petal.Lengthna componente CP1 e das variáveis Sepal. Widthe Sepal.Lengthna segunda. O exemplo teve a finalidade de mostrar a estruturação de uma análise de componentes principal por meio de álgebra matricial e de um conjunto de dados conhecidos. Antes de começarmos a trabalhar e a conhecer o comando princomp(), que é a função do R Base que executa ajustes de componentes principais, sugiro que refaçam todos esses passos com a matriz de correlações, pois, poderá apresentar resultados diferentes quanto à interpretabilidades de suas componentes. Resslta-se que ao utilizarmos as componentes na matriz de covariância de iris, percebemos a dominância da variável Petal.Lengthna análise. Nessa tarefa, foquem em entender cada operação efetuada e o significado dos resultados obtidos.