Métodos Estatísticos Multivariados

Capítulo XVI. Análise de correlação canônica

XVI.1. Combinação linear mais interessante

XVI.1.1 Testando os coeficientes de correlação canônicos

XVI.2. Correlação canônica na prática
XVI.3. Exercícios

Estruturas de dados multivariadas complexas são melhor compreendidas estudando projeções de baixa dimensão. Para um estudo conjunto de dois conjuntos de dados, podemos perguntar que tipo de projeção de baixa dimensão ajuda a encontrar possíveis estruturas conjuntas para as duas amostras. A análise de correlação canônica (CCA) é uma ferramenta padrão de análise estatística multivariada para a descoberta e quantificação de associações entre dois conjuntos de variáveis.

A técnica básica é baseada em projeções. Uma projeção define um índice ou variável multivariada projetada que se correlaciona ao máximo com o índice da outra variável, para cada amostra separadamente. O objetivo do CCA é maximizar a associação, medida por correlação, entre as projeções de baixa dimensão dos dois conjuntos de dados. Os vetores de correlação canônica são encontrados por uma análise de covariância conjunta das duas variáveis. A técnica é aplicada a um exemplo de marketing onde se analisa a associação de um fator preço e outras variáveis, como design, esportividade etc.. Testes são dados sobre como avaliar o significado da associação descoberta.

XVI.1. Combinação linear mais interessante

As associações entre dois conjuntos de variáveis podem ser identificadas e quantificadas por CCA. A técnica foi originalmente desenvolvida por Hotelling (1935) que analisou como a velocidade aritmética e o poder aritmético estão relacionados com a velocidade de leitura e o poder de leitura. Outros exemplos são a relação entre variáveis de política governamental e variáveis de desempenho econômico e a relação entre características do emprego as empresas.

Suponha que recebemos duas variáveis aleatórias \(X\in \mathbb{R}^q\) e \(Y\in\mathbb{R}^p\). A ideia é encontrar um índice que descreva uma possível ligação entre \(X\) e \(Y\). O CCA é baseado em índices lineares, ou seja, combinações lineares \[ a^\top X \qquad \mbox{e} \qquad b^\top Y \]

das variáveis aleatórias.

O CCA procura os vetores \(a\) e \(b\) tais que a relação dos dois índices \(a^\top X\) e \(b^\top Y\) seja quantificada de alguma forma interpretável. Mais precisamente, procura-se as projeções “mais interessantes” \(a\) e \(b\) no sentido de que maximizam a correlação \[ \rho(a,b) = \rho(a^\top X,b^\top Y) \]

entre os dois índices.

Consideremos a correlação \(\rho(a,b)\) entre as duas projeções com mais detalhes. Suponha que \[ \begin{pmatrix} X \\ Y \end{pmatrix} \sim \begin{pmatrix} \begin{pmatrix} \mu \\ \nu \end{pmatrix}, \begin{pmatrix} \Sigma_{XX} & \Sigma_{XY} \\ \Sigma_{YX} & \Sigma_{YY} \end{pmatrix}\end{pmatrix} \]

onde as submatrizes desta estrutura de covariância são dadas por \[ \mbox{Var}(X)=\Sigma_{XX}, \quad \mbox{Var}(Y) = \Sigma_{YY} \] e \[ \mbox{Cov}(X,Y)=\mbox{E}\big( (X-\mu)(Y-\nu)^\top\big) = \Sigma_{XY} = \Sigma_{YX}^\top\cdot \]

Sabemos que \[ \rho(a,b) = \dfrac{a^\top \Sigma_{XY} b}{\big(a^\top\Sigma_{XX}a \big)^{1/2}\big(b^\top\Sigma_{YY} b\big)^{1/2}}\cdot \] Portanto, \(\rho(ca,b)=\rho(a,b)\) para qualquer \(c\in\mathbb{R}^+\).

Dada a invariância da escala, podemos redimensionar as projeções \(a\) e \(b\) e, assim, podemos resolver igualmente \[ \max_{a,b} a^\top \Sigma_{XY} b \] sob as restrições \[ a^\top \Sigma_{XX} a =1 \quad \mbox{e} \quad b^\top \Sigma_{YY} b=1\cdot \]

Para este problema, definamos \[ K = \Sigma_{XX}^{-1/2} \Sigma_{XY} \Sigma_{YY}^{-1/2}\cdot \]

Lembre-se da decomposição em valor singular de \(K\) do Teorema II.2. A matriz \(K\) pode ser decomposta como \[ K = \Gamma \Lambda \Delta^\top \]

com \(\Gamma=(\gamma_1,\cdots,\gamma_k)\), \(\Delta=(\delta_1,\cdots,\delta_k)\) e \(\Lambda=\mbox{diag} (\lambda_1^{1/2},\cdots,\lambda_k^{1/2})\), onde \[ k = \mbox{posto}(K)= \mbox{posto}(\Sigma_{xy})=\mbox{posto}(\Sigma_{YX}) \]

e \(\lambda_1\geq \lambda_2\geq \cdots\geq \lambda_k\) são os autovalores diferentes de zero de \(N_1=KK^\top\) e \(N_2=K^\top K\) e \(\gamma_i\) e \(\delta_j\) são os autovetores padronizados de \(N_1\) e \(N_2\), respectivamente.

Sejam agora, para \(i=1,\cdots,k\) os vetores \[ a_i = \Sigma_{XX}^{-1/2} \gamma_i \qquad \mbox{e} \qquad b_i = \Sigma_{YY}^{-1/2}\delta_i, \]

que são chamados de vetores de correlação canônicos. Usando esses vetores de correlação canônica, definimos as variáveis de correlação canônicas \[ \eta_i = a_i^\top X \qquad \mbox{e} \qquad \varphi_i = b_I^\top Y\cdot \] As quantidades \(\rho_i=\lambda_i^{1/2}\), para \(i = 1,\cdots,k\) são chamados de coeficientes de correlação canônica.

Das propriedades da decomposição de valores singulares temos \[ \mbox{Cov}(\eta_i,\eta_j) = a_i^\top \Sigma_{XX} a_j = \gamma_i^\top \gamma_j =\left\{ \begin{array}{rcl} 1, & \mbox{caso} & i =j, \\ 0, & \mbox{caso} & i\neq j \cdot\end{array}\right. \] O mesmo é verdade para \(\mbox{Cov}(\varphi_i,\varphi_j)\). O seguinte teorema nos diz que os vetores de correlação canônicos são a solução para o problema de maximização.

Teorema XVI.1

Para qualquer \(r\), \(1\leq r\leq k\) o máximo \[ C(r)=\max_{a,b} \, a^\top \Sigma_{XY} b \] restrito a \[ a^\top \Sigma_{XX}a=1, \qquad b^\top \Sigma_{YY}b=1 \qquad \mbox{e} \qquad a_i^\top \Sigma_{XX}a =0, \]

para \(i=1,\cdots,r-1\) é dado por \[ C(r)=\rho_r=\lambda_r^{1/2} \]

e é atingido quando \(a=a_r\) e \(b=b_r\).

Demonstração. A demonstração é em três etapas:

1. Fixando \(a\) e maximizando sobre \(b\), ou seja, resolvendo \[ \max_{b} \, \Big( a^\top \Sigma_{XY}b\Big)^2 = \max_{b} \big(b^\top \Sigma_{YX}a \Big)\Big( a^\top \sigma_{XY}b\Big) \] restrito a \(b^\top\Sigma_{YY}b=1\). Pelo Teorema II.5 o máximo é obtido pela mair autovalorda matriz \[ \Sigma^{-1}_{YY}\Sigma_{YX}aa^\top \Sigma_{XY}\cdot \] Pelo Corolário II.2, o único autovalor diferente de zero é igual \[ a^\top \Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX}a\cdot \]
2. 1. Maximize \(a^\top \Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX}a\) sobre \(a\) restrito às restrições do teorema. Assumindo \(\gamma=\Sigma_{XX}^{1/2} a\) e observando que \(a^\top \Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX}a\) iguala a \[ \gamma^\top \Sigma_{XX}^{-1/2}\Sigma_{XY}\Sigma_{YY}^{-1}\Sigma_{YX}\Sigma {-1/2}_{XX}\gamma = \gamma^\top K^\top K \gamma\cdot \] Assim, resolva o problema equivalente \[ \max_\gamma \, \gamma^\top N_1 \gamma, \] restrito a \(\gamma^\top\gamma=1\), \(\gamma_i^\top\gamma=0\) para \(i=1,\cdots,r-1\).
    Observe que os \(\gamma_i\)’s são os autovetores de \(N_1\) correspondentes aos seus primeiros \(r-1\) maiores autovalores. Assim, como no Teorema XI.3, o máximo em \(\max_\gamma \, \gamma^\top N_1 \gamma\) é obtido igualando \(\gamma\) ao autovetor correspondente ao \(r\)-ésimo maior autovalor, ou seja, \(\gamma=\gamma_r\) ou equivalentemente \(a=a_r\). Isso rende \[ C^2(r)=\gamma_r^\top N_1 \gamma_r =\lambda_r \gamma_r^\top \gamma = \lambda_r\cdot \]

1. Mostrar que o máximo é atingido para \(a=a_r\) e \(b=b_r\). Da decomposição em valores singulares (SVD) de \(K\) concluímos que \(K\delta_r=\rho_r\gamma_r\) e, portanto, \[ a_r^\top \Sigma_{XY}b = \gamma^\top K \delta_r =\rho_r \gamma_r^\top \gamma_r = \rho_r\cdot \]

\(\blacksquare\)

Seja \[ \begin{pmatrix} X \\ Y \end{pmatrix} \sim \begin{pmatrix} \begin{pmatrix} \mu \\ \nu \end{pmatrix}, \begin{pmatrix} \Sigma_{XX} & \Sigma_{XY} \\ \Sigma_{YX} & \Sigma_{YY} \end{pmatrix}\end{pmatrix}\cdot \]

Os vetores de correlação canônica \[ a_1 = \Sigma_{XX}^{-1/2}\gamma_1, \] \[ b_1 = \Sigma_{YY}^{-1/2}\delta_1 \]

maximizam a correlação entre as variáveis canônicas \[ \eta_1 = a_1^\top X, \]

\[ \varphi_1 = b_1^\top Y\cdot \]

A covariância das variáveis canônicas \(\eta\) e \(\varphi\) é dada no próximo teorema.

Teorema XVI.2

Sejam \(\eta_i\) e \(\varphi_i\) sejam as \(i\)-ésimas variáveis de correlação canônicas, \(i=1,\cdots,k\). Definindo \(\eta=(\eta_1,\cdots,\eta_k)\) e \(\varphi=(\varphi_1,\cdots,\varphi_k)\). Então \[ \mbox{Var}\begin{pmatrix} \eta \\ \varphi \end{pmatrix} = \begin{pmatrix} \pmb{I}_k & \Lambda \\ \Lambda & \pmb{I}_k \end{pmatrix}, \] com \(\Lambda=\mbox{diag} (\lambda_1^{1/2},\cdots,\lambda_k^{1/2})\).

Este teorema mostra que os coeficientes de correlação canônica, \(\rho_i=\lambda_i^{1/2}\) são as covariâncias entre as variáveis canônicas \(\eta_i\) e \(\varphi_i\) e que os índices \(\eta_1= a^\top X\) e \(\varphi_1= b^\top Y\) têm a covariância máxima \(\sqrt{\lambda_1}=\rho_1\).

O seguinte teorema mostra que as correlações canônicas são invariantes a transformações lineares das variáveis originais.

Teorema XVI.3

Sejam \(X^* =U^\top X+u\) e \(Y^*= V^\top Y+\nu\) onde \(U\) e \(V\) são matrizes não singulares. Então as correlações canônicas entre \(X^*\) e \(Y^*\) são as mesmas entre \(X\) e \(Y\). Os vetores de correlação canônica de \(X^*\) e \(Y^*\) são dados por \[ a_i^*=U^{-1}a_i, \] \[ b_i^*=V^{-1}b_i\cdot \]

XVI.1.1 Testando os coeficientes de correlação canônicos

A hipótese de que os dois conjuntos de variáveis \(X\) e \(Y\) não são correlacionados pode ser testada, sob suposições de normalidade, com a estatística de razão de verossimilhança de Wilks (Gibbins, 1985): \[ T^{2/n} = \left| \pmb{I}-S_{YY}^{-1} S_{YX}S_{XX}^{-1} S_{XY}\right| = \prod_{1=1}^k \big( 1-\ell_i\big)\cdot \]

Esta estatística infelizmente tem uma distribuição bastante complicada. Bartlett (1939) fornece uma aproximação para \(n\) grande: \[ -\big(n-(p+q+3)/2 \big) \log \left( \prod_{i=1}^k \big( 1-\ell_i\big) \right) \sim \chi^2_{pq}\cdot \]

Um teste da hipótese de que apenas \(s\) dos coeficientes de correlação canônicos são diferentes de zero pode ser baseado (assintoticamente) na estatística \[ -\big(n-(p+q+3)/2 \big) \log \left( \prod_{i=s+1}^k \big( 1-\ell_i\big) \right) \sim \chi^2_{(p-s)(q-s)}\cdot \]

XVI.2. Correlação canônica na prática

Na prática temos que estimar as matrizes de covariância \(\Sigma_{XX}\), \(\Sigma_{XY}\) e \(\Sigma_{YY}\). Vamos aplicar o CCA aos dados das marcas dos carros, CarMarks, Exercício 4, Seção XI.8. No contexto deste conjunto de dados interessa relacionar variáveis de preço com variáveis como esportividade e segurança. Em particular, gostaríamos de investigar a relação entre as duas variáveis de valor e preço não depreciado do carro e todas as outras variáveis.

Exemplo XVI.1

Realizaremos a análise de correlção canônica (CCA) nas matrizes de dados \(X\) e \(Y\) que correspondem ao conjunto de valores \[ \{ \mbox{Price}, \mbox{Value}\} \] e

\[ \{\mbox{Economy}, \mbox{Service}, \mbox{Design}, \mbox{Sporty}, \mbox{Safety}, \mbox{Handling}\}, \] respectivamente.

CarMarks = read.table(file = "http://leg.ufpr.br/~lucambio/MSM/CarMarks.txt", sep = "", header = TRUE)
head(CarMarks)

##      Type    Model Economy Service Value Price Design Sport Safety Handling
## 1    Audi      100     3.9     2.8   2.2   4.2    3.0   3.1    2.4      2.8
## 2     BMW 5 series     4.8     1.6   1.9   5.0    2.0   2.5    1.6      2.8
## 3 Citroen       AX     3.0     3.8   3.8   2.7    4.0   4.4    4.0      2.6
## 4 Ferrari     <NA>     5.3     2.9   2.2   5.9    1.7   1.1    3.3      4.3
## 5    Fiat      Uno     2.1     3.9   4.0   2.6    4.5   4.4    4.4      2.2
## 6    Ford   Fiesta     2.3     3.1   3.4   2.6    3.2   3.3    3.6      2.8

A análise de correlação canônica é usada para identificar e medir as associações entre dois conjuntos de variáveis. A correlação canônica é apropriada nas mesmas situações em que a regressão múltipla seria, mas onde existem múltiplas variáveis de resultado intercorrelacionadas. A análise de correlação canônica determina um conjunto de variáveis canônicas, combinações lineares ortogonais das variáveis dentro de cada conjunto que melhor explicam a variabilidade dentro e entre conjuntos.

library(ggplot2)
library(GGally)
library(CCA)

X = CarMarks[,c(6,5)]
ggpairs(X)

Y = CarMarks[,c(3,4,7,8,9,10)]
ggpairs(Y)

Para vermos as correlações dentro e entre os dois conjuntos de variáveis usamos a função matcor do pacote CCA.

matcor(X, Y)

## $Xcor
##            Price      Value
## Price  1.0000000 -0.8062525
## Value -0.8062525  1.0000000
## 
## $Ycor
##             Economy    Service     Design      Sport     Safety   Handling
## Economy   1.0000000 -0.2669004 -0.6048308 -0.4485906 -0.2699442  0.6060373
## Service  -0.2669004  1.0000000  0.7358654  0.6835027  0.9179337  0.3609505
## Design   -0.6048308  0.7358654  1.0000000  0.8806639  0.7118285 -0.1987946
## Sport    -0.4485906  0.6835027  0.8806639  1.0000000  0.6582605 -0.1995913
## Safety   -0.2699442  0.9179337  0.7118285  0.6582605  1.0000000  0.3322673
## Handling  0.6060373  0.3609505 -0.1987946 -0.1995913  0.3322673  1.0000000
## 
## $XYcor
##               Price      Value    Economy    Service     Design      Sport
## Price     1.0000000 -0.8062525  0.7608951 -0.6820468 -0.8812700 -0.8388292
## Value    -0.8062525  1.0000000 -0.3822370  0.9311772  0.8090157  0.7964358
## Economy   0.7608951 -0.3822370  1.0000000 -0.2669004 -0.6048308 -0.4485906
## Service  -0.6820468  0.9311772 -0.2669004  1.0000000  0.7358654  0.6835027
## Design   -0.8812700  0.8090157 -0.6048308  0.7358654  1.0000000  0.8806639
## Sport    -0.8388292  0.7964358 -0.4485906  0.6835027  0.8806639  1.0000000
## Safety   -0.7045731  0.9019089 -0.2699442  0.9179337  0.7118285  0.6582605
## Handling  0.2926990  0.2478082  0.6060373  0.3609505 -0.1987946 -0.1995913
##              Safety   Handling
## Price    -0.7045731  0.2926990
## Value     0.9019089  0.2478082
## Economy  -0.2699442  0.6060373
## Service   0.9179337  0.3609505
## Design    0.7118285 -0.1987946
## Sport     0.6582605 -0.1995913
## Safety    1.0000000  0.3322673
## Handling  0.3322673  1.0000000

Encontrando e mostrando as correlações canônicas:

CCA = cc(X, Y)
CCA$cor

## [1] 0.9793946 0.9056556

Correlações canônicas brutas:

CCA[3:4]

## $xcoef
##             [,1]      [,2]
## Price -0.3204361 -1.408465
## Value  0.6169295 -1.422128
## 
## $ycoef
##                  [,1]         [,2]
## Economy  -0.434867287 -0.461684146
## Service   0.205117606 -0.702965174
## Design    0.003710477  0.061419728
## Sport     0.467512331 -0.003467896
## Safety    0.221073933  0.295042657
## Handling  0.397165746 -1.006194586

Os coeficientes canônicos brutos são interpretados de maneira análoga à interpretação dos coeficientes de regressão, ou seja, para a variável Economy, um aumento de uma unidade nesta variável leva a uma diminuição de -0.434867287 na primeira variável canônica do conjunto 2 quando todas as outras variáveis são mantidas constantes.

Em seguida, usaremos comput para calcular os carregamentos das variáveis nas dimensões canônicas (variantes). Essas cargas são correlações entre as variáveis e as variáveis canônicas.

# computa carregamentos canônicos
CCA1 = comput(X, Y, CCA)

# exibe carregamentos canônicos
CCA1[3:4]

## $corr.X.xscores
##             [,1]       [,2]
## Price -0.9173969 -0.3979736
## Value  0.9750834 -0.2218386
## 
## $corr.Y.xscores
##                 [,1]       [,2]
## Economy  -0.54248074 -0.6614131
## Service   0.88220655 -0.3198371
## Design    0.87473216  0.1979866
## Sport     0.85035387  0.1475404
## Safety    0.87096398 -0.2373185
## Handling  0.05694861 -0.8667496

As correlações acima são entre variáveis observadas e variáveis canônicas que são conhecidas como cargas canônicas. Essas variáveis canônicas são na verdade um tipo de variável latente.

Em geral, o número de dimensões canônicas é igual ao número de variáveis no conjunto menor; no entanto, o número de dimensões significativas pode ser ainda menor.

As dimensões canônicas, também conhecidas como variáveis canônicas, são variáveis latentes análogas aos fatores obtidos na análise fatorial. Para este modelo em particular, existem três dimensões canônicas, das quais apenas as duas primeiras são estatisticamente significativas.

Para os testes estatísticos usamos o pacote R CCP.

# testa as dimensões canônicas
rho = CCA$cor
rho

## [1] 0.9793946 0.9056556

# Definindo o número de observações, número de variáveis no primeiro conjunto e número de 
# variáveis no segundo conjunto
n = dim(X)[1]
p = length(X)
q = length(Y)

library(CCP)
# Calcula o p-valor utilizando as aproximações F de diferentes testes estatísticos:
p.asym(rho, n, p, q, tstat = "Wilks")

## Wilks' Lambda, using F-approximation (Rao's F):
##                 stat   approx df1 df2      p.value
## 1 to 2:  0.007332857 28.47430  12  32 1.117995e-13
## 2 to 2:  0.179787957 15.51117   5  17 8.143198e-06

p.asym(rho, n, p, q, tstat = "Hotelling")

##  Hotelling-Lawley Trace, using F-approximation:
##               stat   approx df1 df2      p.value
## 1 to 2:  28.080239 35.10030  12  30 2.464695e-14
## 2 to 2:   4.562108 15.51117   5  34 5.954839e-08

p.asym(rho, n, p, q, tstat = "Pillai")

##  Pillai-Bartlett Trace, using F-approximation:
##              stat    approx df1 df2      p.value
## 1 to 2:  1.779426 22.857202  12  34 8.031353e-13
## 2 to 2:  0.820212  5.283671   5  38 8.828018e-04

p.asym(rho, n, p, q, tstat = "Roy")

##  Roy's Largest Root, using F-approximation:
##               stat   approx df1 df2      p.value
## 1 to 1:  0.9592139 66.63471   6  17 7.160728e-11
## 
##  F statistic for Roy's Greatest Root is an upper bound.

A variável canônica \(\widehat{\eta}_1\) pode ser interpretada como um índice de Price e Value. A variável canônica \(\widehat{\varphi}_1\) é formada principalmente pelas variáveis Service, Design, Sport e Safety.

XVI.3. Exercícios

Mostre que os autovalores de \(KK^\top\) e \(K^\top K\) são idênticos.

Faça o CCA para os seguintes subconjuntos de variáveis: \(X\) correspondente a fpriceg e Y correspondente a feconomy, easy handlingg a partir dos dados de marcas de carros (Tabela 22.7).

Calcule as primeiras variáveis canônicas para o Exemplo 16.1. Interprete os coeficientes.

Use a decomposição em valores proprios da matriz \(K\) para mostrar que as variáveis canônicas \(\eta_1\) e \(\eta_2\) não são correlacionadas.

Verifique se o número de autovalores diferentes de zero da matriz \(K\) é igual a \(\mbox{posto}(\Sigma_{XY})\).

Expresse a decomposição de valores singulares ou proprios das matrizes \(K\) e \(K^\top\) usando autovalores e autovetores das matrizes \(K^\top K\) e \(KK^\top\).

Qual será o resultado do CCA para \(Y =X\)?

Quais serão os resultados do CCA para \(Y=2X\) e para \(Y=-X\)?

Que resultados você espera se executar CCA para \(X\) e \(Y\) de modo que \(\Sigma_{XY}=0\)? E se \(\Sigma_{XY}=\pmb{I}_p\)?