20 Análise de Correlação Canônica com Pinguins

Neste exemplo, aplicamos a Análise de Correlação Canônica (ACC) para investigar a relação entre medidas do bico e medidas corporais dos pinguins do arquipélago Palmer. Utilizamos Python e o dataset palmerpenguins, que contém medidas de três espécies (Adelie, Chinstrap, Gentoo).

Conjuntos de variáveis:

Grupo 1 - Bico (\(\mathbf{X}^{(1)}\)): bill_length_mm (comprimento, na direção do bico) e bill_depth_mm (profundidade/altura, medida vertical do bico)
Grupo 2 - Corpo (\(\mathbf{X}^{(2)}\)): flipper_length_mm (nadadeira) e body_mass_g (massa)

Objetivo: Quantificar a associação entre a forma do bico e o tamanho corporal, e interpretar essa relação.

20.1 Preparação e Análise Exploratória

Código

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from scipy import linalg
from scipy.stats import chi2 as chi2_dist

sns.set_theme(style="whitegrid")

# Carregar e limpar dados
penguins = sns.load_dataset("penguins")
penguins_clean = penguins.dropna()

# Definir os dois grupos
X = penguins_clean[['bill_length_mm', 'bill_depth_mm']]
Y = penguins_clean[['flipper_length_mm', 'body_mass_g']]

# Padronizar (trabalharemos com a matriz de correlação)
X_std = (X - X.mean()) / X.std()
Y_std = (Y - Y.mean()) / Y.std()

print(f"Amostra: n = {len(penguins_clean)}, p = {X.shape[1]}, q = {Y.shape[1]}")

Amostra: n = 333, p = 2, q = 2

Antes da ACC, visualizamos a matriz de correlação para verificar as relações entre e dentro dos grupos:

Código

combined = pd.concat([X_std, Y_std], axis=1)
corr_matrix = combined.corr()

plt.figure(figsize=(8, 6))
sns.heatmap(corr_matrix, annot=True, cmap='coolwarm', vmin=-1, vmax=1, fmt=".2f", 
            cbar_kws={'label': 'Correlação'})
plt.title("Correlações Intra e Inter-Grupos")
plt.tight_layout()
plt.show()

Figura 20.1: Matriz de correlação entre as variáveis dos dois grupos.

Observamos correlações cruzadas substanciais: flipper_length_mm correlaciona positivamente com bill_length_mm (0.66) e negativamente com bill_depth_mm (-0.58). A ACC sintetizará essas relações.

20.2 Cálculo da ACC via SVD

Seguimos a abordagem via SVD apresentada em Capítulo 12. Particionamos a matriz de correlação \(\mathbf{R}\) e aplicamos SVD à matriz de coerência \(\mathbf{K} = \mathbf{R}_{11}^{-1/2} \mathbf{R}_{12} \mathbf{R}_{22}^{-1/2}\):

Código

R = corr_matrix.values
p, q = X.shape[1], Y.shape[1]

# Partições da matriz de correlação
R11 = R[:p, :p]
R22 = R[p:, p:]
R12 = R[:p, p:]
R21 = R[p:, :p]

# Branqueamento: inversas das raízes quadradas
R11_inv_sqrt = linalg.fractional_matrix_power(R11, -0.5)
R22_inv_sqrt = linalg.fractional_matrix_power(R22, -0.5)

# Matriz de coerência K
K = R11_inv_sqrt @ R12 @ R22_inv_sqrt

# SVD de K
U_k, canonical_corrs, Vh_k = linalg.svd(K)

# Coeficientes canônicos (revertendo o branqueamento)
A = R11_inv_sqrt @ U_k
B = R22_inv_sqrt @ Vh_k.T

# Ajuste de sinal para interpretação intuitiva:
# SVD pode retornar vetores com sinal arbitrário. Para facilitar a interpretação,
# invertemos os sinais de modo que a primeira variável do corpo (flipper_length) 
# tenha carga positiva na primeira dimensão canônica.
if (R22 @ B[:, 0])[0] < 0:
    A[:, 0] = -A[:, 0]
    B[:, 0] = -B[:, 0]

print("Correlações Canônicas (\rho^*):")
for i, rho in enumerate(canonical_corrs, 1):
    print(f"  rho_{i}* = {rho:.4f}")

Correlações Canônicas (
ho^*):
  rho_1* = 0.7876
  rho_2* = 0.0864

As duas correlações canônicas obtidas (\(\rho_1^* = 0.7876\) e \(\rho_2^* = 0.0864\)) revelam que apenas a primeira dimensão é substancial (> 0.3), conforme discutido em Capítulo 12. A segunda é desprezível, sugerindo que a relação entre os dois grupos se concentra essencialmente em uma única dimensão latente.

Nota sobre Sinais das Variáveis Canônicas

Os sinais dos vetores canônicos retornados pela SVD são arbitrários - multiplicar ambos \(U_k\) e \(V_k\) por -1 resulta na mesma correlação canônica. Para facilitar a interpretação, ajustamos os sinais de modo que valores positivos de \(V_1\) correspondam a pinguins maiores.

20.3 Interpretação via Cargas Canônicas

Conforme a teoria, os coeficientes \(\mathbf{a}_k\) e \(\mathbf{b}_k\) são difíceis de interpretar devido ao branqueamento. Calculamos as cargas canônicas (correlações entre variáveis originais e canônicas) e cargas cruzadas (correlações entre variáveis de um grupo e canônicas do outro):

Código

# Cargas canônicas (l_X,k e l_Y,k)
loadings_X = R11 @ A  # Corr(X^(1), U)
loadings_Y = R22 @ B  # Corr(X^(2), V)

# Cargas cruzadas (l_X->Y,k e l_Y->X,k)
cross_loadings_X = R12 @ B  # Corr(X^(1), V)
cross_loadings_Y = R21 @ A  # Corr(X^(2), U)

# Criar DataFrames para visualização
df_loadings_X = pd.DataFrame(loadings_X, index=X.columns, columns=['U1', 'U2'])
df_loadings_Y = pd.DataFrame(loadings_Y, index=Y.columns, columns=['V1', 'V2'])

print("Cargas Canônicas - Primeira Dimensão:")
print("\nGrupo Bico (correlação com U_1):")
print(df_loadings_X['U1'].round(3))
print("\nGrupo Corpo (correlação com V_1):")
print(df_loadings_Y['V1'].round(3))

Cargas Canônicas - Primeira Dimensão:

Grupo Bico (correlação com U_1):
bill_length_mm    0.828
bill_depth_mm    -0.735
Name: U1, dtype: float64

Grupo Corpo (correlação com V_1):
flipper_length_mm    1.000
body_mass_g          0.866
Name: V1, dtype: float64

Visualização das cargas:

Código

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

# Grupo Bico
sns.barplot(x=df_loadings_X['U1'], y=df_loadings_X.index, ax=axes[0], 
            color='steelblue', edgecolor='black')
axes[0].set_title('Cargas Canônicas: U_1 (Bico)', fontsize=12, fontweight='bold')
axes[0].set_xlabel('Correlação com U_1')
axes[0].set_ylabel('')
axes[0].set_xlim(-1, 1)
axes[0].axvline(0, color='black', linewidth=0.8)

# Grupo Corpo
sns.barplot(x=df_loadings_Y['V1'], y=df_loadings_Y.index, ax=axes[1], 
            color='coral', edgecolor='black')
axes[1].set_title('Cargas Canônicas: V_1 (Corpo)', fontsize=12, fontweight='bold')
axes[1].set_xlabel('Correlação com V_1')
axes[1].set_ylabel('')
axes[1].set_xlim(-0.5, 1.1)
axes[1].axvline(0, color='black', linewidth=0.8)

plt.tight_layout()
plt.show()

Figura 20.2: Cargas canônicas da primeira dimensão. Valores indicam correlação entre variáveis originais e variáveis canônicas.

Interpretação:

Observando as cargas canônicas:

Dimensão Canônica do Bico (\(U_1\)): bill_length_mm (comprimento do bico) tem carga positiva elevada (+0.83), enquanto bill_depth_mm (altura/espessura do bico) tem carga negativa (-0.74). Isso define um contraste morfológico:
- \(U_1\) positivo -> bicos longos e finos (alto comprimento, baixa profundidade)
- \(U_1\) negativo -> bicos curtos e grossos (baixo comprimento, alta profundidade)
Dimensão Canônica do Corpo (\(V_1\)): Ambas as variáveis têm cargas positivas elevadas (flipper_length_mm: +1.00, body_mass_g: +0.87), representando o tamanho geral.
- \(V_1\) positivo -> pinguins maiores
- \(V_1\) negativo -> pinguins menores

A correlação canônica de 0.79 indica forte associação positiva: pinguins com bicos longos e finos (\(U_1\) alto) tendem a ser maiores (\(V_1\) alto), enquanto pinguins com bicos curtos e grossos (\(U_1\) baixo) tendem a ser menores (\(V_1\) baixo).

Índice de redundância: Conforme discutido em Capítulo 12, o índice de redundância mede quanto da variância de um grupo é explicada pela variável canônica do outro grupo:

Código

# Variância extraída (média dos quadrados das cargas cruzadas)
var_ext_X = np.mean(cross_loadings_X[:, 0]**2)
var_ext_Y = np.mean(cross_loadings_Y[:, 0]**2)

# Redundância = Variância Extraída × \rho^2
red_X = var_ext_X * (canonical_corrs[0]**2)
red_Y = var_ext_Y * (canonical_corrs[0]**2)

print(f"Grupo Bico: Var. Extraída = {var_ext_X:.3f}, Redundância = {red_X:.3f}")
print(f"Grupo Corpo: Var. Extraída = {var_ext_Y:.3f}, Redundância = {red_Y:.3f}")

Grupo Bico: Var. Extraída = 0.380, Redundância = 0.236
Grupo Corpo: Var. Extraída = 0.543, Redundância = 0.337

Interpretação da redundância:

Variância Extraída: Indica quanto da variância de um grupo é capturada pela correlação com a variável canônica do outro grupo. Por exemplo, 54.3% da variância das medidas corporais (\(V_1\)) correlaciona com as medidas do bico (\(U_1\)).
Redundância: Ajusta a variância extraída pela correlação canônica ao quadrado. Aqui, a dimensão canônica do bico (\(U_1\)) explica 33.7% da variância total das medidas corporais (0.543 × 0.79² = 0.337). Reciprocamente, \(V_1\) explica apenas 23.6% da variância das medidas do bico, indicando que o tamanho corporal é mais previsível a partir da forma do bico do que o inverso.

20.4 Visualização no Espaço Canônico

Calculamos os scores das variáveis canônicas para cada observação e visualizamos a distribuição das espécies:

Código

# Scores canônicos
U = X_std @ A
V = Y_std @ B

scores_df = pd.DataFrame({
    'U1': U.iloc[:, 0],
    'V1': V.iloc[:, 0],
    'Species': penguins_clean['species']
})

plt.figure(figsize=(8, 5))
sns.scatterplot(data=scores_df, x='U1', y='V1', hue='Species', style='Species', 
                s=100, alpha=0.7, palette='Set2')

# Linha representando a correlação canônica
x_vals = np.array([scores_df['U1'].min(), scores_df['U1'].max()])
plt.plot(x_vals, x_vals * canonical_corrs[0], color='gray', linestyle='--', 
         linewidth=2, label=f'rho_1* = {canonical_corrs[0]:.2f}')

plt.xlabel('U_1 (Bico: Curto/Grosso <-> Longo/Fino)', fontsize=11)
plt.ylabel('V_1 (Corpo: Menor <-> Maior)', fontsize=11)
plt.title('Espaço Canônico: Forma do Bico vs. Tamanho Corporal', fontsize=12, fontweight='bold')
plt.legend(title='Espécie', fontsize=10)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

Figura 20.3: Distribuição das espécies no espaço das variáveis canônicas. A correlação canônica de 0.79 é evidente na linearidade, e a separação das espécies revela diferenças morfológicas.

O gráfico confirma:

Correlação linear forte entre \(U_1\) e \(V_1\) (linha tracejada)
Gentoo (maiores, bicos longos/finos) no quadrante superior direito
Adelie (menores, bicos curtos/grossos) no quadrante inferior esquerdo
Chinstrap em posição intermediária

20.5 Teste de Significância

Conforme Capítulo 12, testamos \(H_0: \mathbf{\Sigma}_{12} = \mathbf{0}\) usando a estatística Lambda de Wilks e a aproximação qui-quadrado de Bartlett:

Código

n = len(penguins_clean)

# Lambda de Wilks
wilks_lambda = np.prod(1 - canonical_corrs**2)

# Estatística qui-quadrado de Bartlett
chi2_stat = -(n - 1 - (p + q + 1) / 2) * np.log(wilks_lambda)
df = p * q
p_value = 1 - chi2_dist.cdf(chi2_stat, df)

print(f"Lambda de Wilks (Lambda): {wilks_lambda:.4f}")
print(f"Qui-quadrado: chi^2 = {chi2_stat:.2f}, gl = {df}")
print(f"P-valor: {p_value:.2e}")
print(f"\nConclusão: {'Rejeita H_0' if p_value < 0.05 else 'Não rejeita H_0'} (alpha = 0.05)")

Lambda de Wilks (Lambda): 0.3768
Qui-quadrado: chi^2 = 321.60, gl = 4
P-valor: 0.00e+00

Conclusão: Rejeita H_0 (alpha = 0.05)

O p-valor extremamente baixo (\(p < 0.001\)) rejeita a hipótese de independência. A associação entre forma do bico e tamanho corporal é estatisticamente significativa.

Síntese: A ACC revelou uma dimensão canônica dominante (\(\rho_1^* = 0.79\)) que relaciona a morfologia do bico (contraste comprimento/profundidade) ao tamanho corporal geral. Essa estrutura reflete diferenças ecológicas entre as espécies de pinguins, com Gentoo exibindo um fenótipo distinto (maior porte, bicos alongados) comparado a Adelie.