- Warning 1: This post is written in Portuguese.
- Warning 2: Most of what I have here I got from other posts/examples online. All links are provided and I encourage you to take a look at them.
A Análise de Componentes Principais (em inglês PCA) é o nome comum dado à técnica que usa princípios de álgebra linear para transformar variáveis, possivelmente correlacionadas, em um número menor de variáveis chamadas de Componentes Principais (novamente em inglês PC).
A PCA é usada em diversas aplicações, desde a compressão de dados (MP3, JPG) até remoção de ruídos, passando pela análises de grande quantidade de dados.
PCA do ponto de vista da Geometria Espacial
Em termos gerais a PCA busca reduzir o número de dimensões de um set de dados. Projetando os dados em um novo plano. Usando essa nova projeção os dados originais, que podem envolver diversas variáveis, podem ser interpretados utilizando menos "dimensões."
No set de dados reduzido podemos observar com mais clareza tendências, padrões e/ou outliers. Mas vale lembrar que a regra: "Se não está nos dados brutos não existe!" é sempre válida. A PCA fornece apenas mais clareza aos padrões que já estão lá.
Utilizaremos dois vetores 1D conhecidos $(x, y)$ para entender essa nova projeção que a PCA faz. Note que os vetores possuem uma clara dependência linear.
O primeiro passo é remover as médias dos dados.
(Dados e exemplo: http://bekoc.blogspot.com.br/2013/12/implementing-principle-component.html)
import numpy as np
from pandas import DataFrame
x = np.array([2.5, 0.5, 2.2, 1.9, 3.1, 2.3, 2.0, 1.0, 1.5, 1.1])
y = np.array([2.4, 0.7, 2.9, 2.2, 3.0, 2.7, 1.6, 1.1, 1.6, 0.9])
Z = np.c_[x - x.mean(), y - y.mean()]
df = DataFrame(Z, columns=['x', 'y'])
df.index.name = 'Medidas'
df.T
O segundo passo é a normalização dos dados. A mais costumeira é, após remover a média, dividir os dados pelo desvio padrão. Esse procedimento também é chamado de z-score.
import matplotlib.pyplot as plt
line = dict(linewidth=1, linestyle='--', color='k')
marker = dict(linestyle='none', marker='o', markersize=7, color='blue', alpha=0.5)
fig, ax = plt.subplots(figsize=(5, 5))
ax.plot(x, y, **marker)
ax.axhline(**line)
ax.axvline(**line)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title(u'Dados Originais (com a média)')
_ = ax.axis([-1, 4, -1, 4])
O terceiro passo é calcular a covariância dos dados normalizados.
cov = np.cov(Z.T)
cov
A "redução" dos dados pode ser feita por Singular Value Decomposition (SVD) ou Empirical Orthogonal Functions (EOF). Estamos buscando planos ortogonais entre as variáveis que maximize a variança.
Vamos utilizar a EOF, ou seja, calculando os Autovalores e Autovetores
da matriz de covariança. Mas colocarei aqui o correspondente usando a SVD como exercício:
U, S, V = np.linalg.svd(cov, full_matrices=True, compute_uv=True)
S # Similar aos Autovalores
V # Similar aos Autovetores
PCs = np.dot(V, Z.T) # Crie um gráfico que se compare com os abaixo.
eigenvalues, eigenvectors = np.linalg.eig(cov)
eigenvalues
eigenvectors
As componentes principais são os Autovetores. É comum normalizar os Autovetores para facilitar as comparações entre eles.
pc = eigenvectors[1] / eigenvectors[0]
pc
Vamos plotar os vetores que definem esse novo plano ortogonal saindo do ponto $(0, 0)$ até $x$ máximo (Autovalor) e $y$ máximo (componente principal normalizada e dimensionada para o autovalor.)
fig, ax = plt.subplots(figsize=(5, 5))
ax.plot(Z[:, 0], Z[:, 1], **marker)
ax.set_xlim([-2, 2])
ax.set_ylim([-2, 2])
ax.axhline(**line)
ax.axvline(**line)
ax.set_xlabel('x')
ax.set_ylabel('y')
arrowprops = dict(width=0.01, head_width=0.05, alpha = 0.5,
length_includes_head=False)
a1 = ax.arrow(0, 0, eigenvalues[0], pc[0] * eigenvalues[0],
color='k', **arrowprops)
a2 = ax.arrow(0, 0, eigenvalues[1], pc[1] * eigenvalues[1],
color='k', **arrowprops)
a3 = ax.arrow(0, 0, -eigenvalues[0], -pc[0] * eigenvalues[0],
color='r', **arrowprops)
a4 = ax.arrow(0, 0, -eigenvalues[1], -pc[1] * eigenvalues[1],
color='r', **arrowprops)
É fácil de ver que esse novo plano se orienta onde os dados variam. Bom, geometricamente e isso. Agora vamos testar um pacote pronto para a análise de PCA. Vou comparar o que acabamos de fazer com os resultados do módulo scikit-learn.
from sklearn.decomposition import PCA
pca = PCA(n_components=2, copy=True)
X = pca.fit_transform(Z)
fig, ax = plt.subplots()
ax.plot(X[:, 0], X[:, 1], **marker)
ax.set_xlim([-2, 2])
ax.set_ylim([-2, 2])
ax.axhline(**line)
ax.axvline(**line)
_ = ax.set_title("PCA")
Essa seria a cara dos dados "rodados" para o novo plano calculado pela análise. Para entender melhor o que isso quer dizer de uma olhada em:
Note que os dados em si são os mesmos. Nos apenas rodamos para um novo set de eixos (principais). Estamos literalmente olhado para os dados sob um novo ângulo ;-). Esse novo "ângulo" é mais intuitivo para tirarmos conclusões sobre os dados.
Um exemplo de aplicação direta em oceanografia física é rodar as componentes da velocidade $u$ e $v$ de um fundeio em rio ou na plataforma continental para os seus eixos principais. Intuitivamente podemos concluir que tais eixos, em geral, corresponderão aos "ao longo" e "perpendicular" à margem do rio ou dá isóbata do fundeio costeiro. Não é mais intuitivo que eixos geográficos Norte e Sul?
(Aqui tem uma função faz exatamente isso para séries temporais de velocidade.)
Vamos comprar os resultados entre sklearn e o que fizemos anteriormente:
print(eigenvectors)
print('')
print(pca.components_)
Mesmo resultado! Mas note que as componentes estão organizadas
em ordem crescente dos autovalores no módulo sklearn. O que
mais temos no objeto pca:
pca.explained_variance_ratio_ # Quanto da variância é explicada por cada componente.
Para finalizar essa parte deixo
aqui
um exemplo 3D que compara outros módulos além do sklearn.
Entendendo a PCA do ponto de visto estatístico
Uma PCA entre duas variáveis, como fizemos no exemplo acima, não é muito diferente de uma correlação simples entre elas. A PCA começa a ser útil quando lidamos com diversas variáveis. Para entender essa afirmação vamos calcular a correlação de Pearson para os nossos dados.
from scipy.stats.stats import pearsonr
r, p = pearsonr(x-x.mean(), y-y.mean())
print('Pearson r: {}\nPC 1: {}:'.format(r, pc[0]))
Bem próximo da variância explicada pela primeira PC. Já que estamos tentado relacionar esse resultado com correlação linear simples, como seria a relação entre a primeira Componente Principal e um ajuste linear?
import statsmodels.api as sm
def OLSreg(y, Xmat):
return sm.OLS(y, sm.add_constant(Xmat, prepend=True)).fit()
scale = lambda x: (x - x.mean()) / x.std()
ols_fit = OLSreg(Z[:, 1], Z[:, 0])
fig, ax = plt.subplots(figsize=(5, 5))
ax.plot(Z[:, 0], Z[:, 1], **marker)
ax.set_xlim([-2, 2])
ax.set_ylim([-2, 2])
ax.axhline(**line)
ax.axvline(**line)
ax.set_xlabel('x')
ax.set_xlabel('y')
ax.plot(Z[:, 0], ols_fit.fittedvalues, 'r', alpha=0.5)
ax.text(-1, 1, r'R$^2$ = %4.3f' % round(ols_fit.rsquared, 3))
a1 = ax.arrow(0, 0, eigenvalues[1], pc[1] * eigenvalues[1],
**arrowprops)
a2 = ax.arrow(0, 0, -eigenvalues[1], -pc[1] * eigenvalues[1],
**arrowprops)
Há um ângulo pequeno entre as retas formadas pela correlação linear (vermelho) e a reta da primeira Componente Principal (azul). Eu gosto de pensar que a correlação linear está "contaminada" pela variância da segunda Componente Principal.
PCA com dados reais
Para realmente entendermos o que é a PCA temos que usar dados reais. Escolhi seguir o exemplo disponível nesse PDF que usa dados de consumo de comida em gramas, por pessoa, por semana, no Reino Unido.
%%file data.csv
Food,England,Wales,Scotland,N Ireland
Cheese,105,103,103,66
Carcass meat,245,227,242,267
Other meat,685,803,750,586
Fish,147,160,122,93
Fats and oils,193,235,184,209
Sugars,156,175,147,139
Fresh potatoes,720,874,566,1033
Fresh Veg,253,265,171,143
Other Veg,488,570,418,355
Processed potatoes,198,203,220,187
Processed Veg,360,365,337,334
Fresh fruit,1102,1137,957,674
Cereals,1472,1582,1462,1494
Beverages,57,73,53,47
Soft drinks,1374,1256,1572,1506
Alcoholic drinks,375,475,458,135
Confectionery,54,64,62,41
from pandas import read_csv
df = read_csv('data.csv', index_col='Food')
df
Sem realizar ajustes lineares para cada par de dados fica quase impossível
visualizar algum padrão ou tendência nos dados acima. Esse é um caso onda a
PCA pode ajudar. Novamente vamos usar o módulo sklearn. Vale lembrar que
esse módulo usa a SVD por trás da cortina, mas nos poupa do trabalho de lidar
as saídas U, S e V já calculando as propriedades que queremos
diretamente.
Começaremos definindo uma função para normalizar os dados.
def z_score(x):
"""Remove a média e normaliza os pelo desvio padrão"""
return (x - x.mean()) / x.std()
from sklearn.decomposition import PCA
pca = PCA(n_components=None)
pca.fit(df.apply(z_score).T)
Para facilitar vamos criar um um novo DataFrame com os resultados do objeto
pca.
from pandas import DataFrame
loadings = DataFrame(pca.components_.T)
loadings.index = ['PC %s' % pc for pc in loadings.index + 1]
loadings.columns = ['TS %s' % pc for pc in loadings.columns + 1]
loadings
Dependendo dos dados fica mais fácil interpretar as componentes principais quando essas são projetadas nos dados originais. Isso faz com que os números tenham a mesma unidades dos dados originais.
Atenção! Alguns papers não projetam. Fique atento(a) para saber o que está lendo nos gráficos que encontrar por aí.
PCs = np.dot(loadings.values.T, df)
Vamos ignorar a segunda dimensão (PC 2 = 0) e plotar apenas a primeira.
marker = dict(linestyle='none', marker='o', markersize=7, color='blue', alpha=0.5)
fig, ax = plt.subplots(figsize=(7, 2.75))
ax.plot(PCs[0], np.zeros_like(PCs[0]),
label="Scores", **marker)
[ax.text(x, y, t) for x, y, t in zip(PCs[0], loadings.values[1, :], df.columns)]
ax.set_xlabel("PC1")
_ = ax.set_ylim(-1, 1)
marker = dict(linestyle='none', marker='o', markersize=7, color='blue', alpha=0.5)
Facilmente vemos que os habitantes da Irlanda do Norte têm uma dieta diferente da Inglaterra, Escócia e País de Gales. Esses três últimos parecem se agrupar em um grupo de dieta similar. (O que seria uma boa hipótese inicial, já que a Irlanda do Norte está geograficamente e culturalmente mais afastada do resto do Reino Unido. A propósito, muitas das hipóteses que vemos em papers são definidas assim, após as análises ;-)
Esse resultado não é tão evidente quando usamos uma correlação simples. Note que a Irlanda do Norte realmente se correlaciona menos com o resto, mas não tão fortemente como vemos na PC 1 acima. Isso ocorre porque temos a "contaminação" da PC 2 nessa correlação e, nesse caso, a PC 2 contém uma boa porção da variância como veremos mais a frente.
ax = seaborn.corrplot(df, annot=True, diag_names=False)
Agora vamos plotar a primeira e segunda PCs juntas. Esse tipo de gráfico é chamado de Score plot.
fig, ax = plt.subplots(figsize=(7, 2.75))
ax.plot(PCs[0], PCs[1], label="Scores", **marker)
ax.set_xlabel("PC1")
ax.set_ylabel("PC2")
text = [ax.text(x, y, t) for x, y, t in
zip(PCs[0], PCs[1]+0.5, df.columns)]
Note que na "segunda" dimensão estão as diferenças entre os três países que agrupamos no gráfico anterior.
Vamos criar mais alguns gráficos úteis para ajudar na interpretação dos dados. Começaremos com a pergunta: o quanto da variância dos dados é explicada por cada componente?
(Isso nos informa se estamos OK usando 1, 2 ou mais dimensões. Em geral se busca componentes o suficiente para explicar entre 70-80% dos dados.)
perc = pca.explained_variance_ratio_ * 100
perc = DataFrame(perc, columns=['Percentage explained ratio'], index=['PC %s' % pc for pc in np.arange(len(perc)) + 1])
ax = perc.plot(kind='bar')
As componentes 1 e 2 juntas explicam mais de 96% da variância dos dados. Saímos de uma espaço de 17x4 para um espaço de 2x4! (Se isso fosse compressão de uma foto para um arquivo JPEG nos teríamos um arquivo bem reduzido.)
Vamos olhar também a "série" para a primeira componente.
TS1 = loadings['TS 1']
TS1.index = df.index
ax = TS1.plot(kind='bar')
Note que Fresh Potatoes é o inverso de Fresh Fruit, certo? Checando os dados originais vemos que a Irlanda do Norte consome mais Fresh Potatoes e menos Fresh Fruit que o resto do Reino Unido. A diferença entra Irlanda do Norte e o resto não foi o resultado mais notável da primeira componente? As coisas começam a se encaixar!
Outro gráfico comum para explorar os resultados é o Loadings plot ou seja, a influência de cada variável original nas componentes principais. Note que Fresh Fruit, Alcoholic drinks, Soft drinks e Fresh potatoes se destacam da aglomeração central.
marker = dict(linestyle='none', marker='o', markersize=7, color='blue', alpha=0.5)
fig, ax = plt.subplots(figsize=(7, 7))
ax.plot(loadings.icol(0), loadings.icol(1), label="Loadings", **marker)
ax.set_xlabel("non-projected PC1")
ax.set_ylabel("non-projected PC2")
ax.axis([-1, 1, -1, 1])
text = [ax.text(x, y, t) for
x, y, t in zip(loadings.icol(0), loadings.icol(1), df.index)]
Fica para exercício do leitor olhar novamente a tabela original dos dados e procurar algum padrão entre os países e as comidas no Loadings Plot.
Depois disso podemos retornar aos dados originais e, de posse desse conhecimento, começar e discutir os dados.
HTML(html)
This post was written as an IPython notebook. It is available for download or as a static html.
python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.