Análise de Correlação Canônica

Anderson Rodrigo da Silva. (2016). Métodos de Análise Multivariada em R. FEALQ. Piracicaba.

Aplicações em:

• Num estudo de avaliação nutricional de plantas, pode ser de interesse saber se há alguma associação entre os nutrientes (variáveis) encontrados nas raízes (grupo \(1\)) e na parte aérea da planta (grupo \(2\));
• Pode ser de interesse verificar se variáveis químicas (grupo \(1\)) de determinado tipo de fruto são correlacionadas com variáveis morfométricas (grupo \(2\));
• Um biólogo pode estar interessado em saber se variáveis biométricas (grupo \(1\)) de determinada espécie de peixe estão relacionadas a variáveis de qualidade da água (grupo \(2\));
• Características sensoriais (grupo \(1\)) de determinado alimento após certo tempo de exposição ao ambiente, estão relacionadas a variáveis ambientais (grupo \(2\));
• Variáveis de consumo/durabilidade (grupo \(1\)) de determinado veículo estão relacionadas com variáveis socioeconômicas (grupo \(2\)) dos seus donos.

Exemplo 12.1

Referência: Daniel Ferreira Furtado 3a. Ed. Editora UFLA

A partir da matriz de covariâncias amostral dos solos de pastagem da Amazônia, considerando as variáveis Ca e Mg (Cálcio e Magnésio) no primeiro grupo e as variáveis SB e t (Saturação de Bases e Capacidade de troca de cátions) no segundo. A correlação entre os dois grupos de variáveis é utilizada para descrever as relações de causa e efeito em relação à acidez dos solos.

• Estimar as correlações canônicas; e
• Determinar as variáveis canônicas

A amostra de \(n=13\) pontos é dada por:

bd<-matrix(c(
  2.3 ,1.7 ,4.1 ,7.5,
  2.5 ,2.5 ,5.1 ,7.4,
  1.8 ,2.1 ,4.1 ,6.4,
  3.4 ,2.5 ,6.1 ,8.4,
  1.8 ,1.1 ,3.0 ,6.5,
  3.7 ,1.4 ,5.2 ,8.5,
  1.4 ,0.7 ,2.2 ,4.7,
  1.5 ,0.6 ,2.2 ,4.0,
  2.8 ,2.2 ,5.1 ,9.0,
  1.4 ,0.8 ,2.3 ,3.7,
  1.8 ,0.6 ,2.5 ,7.2,
  1.9 ,1.7 ,3.7 ,6.0,
  2.8 ,0.8 ,3.7 ,5.9),nrow=13,byrow = TRUE)

colnames(bd)<-c("Ca", "Mg", "SB", "t");

Variáveis Canônicas e Correlação Canônica Amostrais

Estimação do Parâmetros

Banco de dados: solos de pastagem da Amazônia

bd

       Ca  Mg  SB   t
 [1,] 2.3 1.7 4.1 7.5
 [2,] 2.5 2.5 5.1 7.4
 [3,] 1.8 2.1 4.1 6.4
 [4,] 3.4 2.5 6.1 8.4
 [5,] 1.8 1.1 3.0 6.5
 [6,] 3.7 1.4 5.2 8.5
 [7,] 1.4 0.7 2.2 4.7
 [8,] 1.5 0.6 2.2 4.0
 [9,] 2.8 2.2 5.1 9.0
[10,] 1.4 0.8 2.3 3.7
[11,] 1.8 0.6 2.5 7.2
[12,] 1.9 1.7 3.7 6.0
[13,] 2.8 0.8 3.7 5.9

Tamanho amostral

n<-13

Vetor de Média amostral

colMeans(bd)

      Ca       Mg       SB        t 
2.238462 1.438462 3.792308 6.553846 

Matriz de variância-covariancia estimada

s<-cov(bd);s

          Ca        Mg        SB         t
Ca 0.5692308 0.2875641 0.8628205 1.0085897
Mg 0.2875641 0.5242308 0.8261538 0.8227564
SB 0.8628205 0.8261538 1.7107692 1.8454487
t  1.0085897 0.8227564 1.8454487 2.8393590

Matriz de correlação estimada

r<-cor(bd);r

          Ca        Mg        SB         t
Ca 1.0000000 0.5264162 0.8743404 0.7933416
Mg 0.5264162 1.0000000 0.8723765 0.6743725
SB 0.8743404 0.8723765 1.0000000 0.8373290
t  0.7933416 0.6743725 0.8373290 1.0000000

No. variaveis do \(1^o\) grupo

p<-2 

No. variaveis do \(2^o\) grupo

q<-2 

Número de variáveis canônicas

nvc<-min(p,q);nvc 

[1] 2

Covariância do \(1^o\) grupo: \(\mathbf{s}_{11}\)

s11<-s[1:2,1:2];s11

          Ca        Mg
Ca 0.5692308 0.2875641
Mg 0.2875641 0.5242308

Covariância do \(2^o\) grupo: \(\mathbf{s}_{22}\)

s22<-s[3:4,3:4];s22

         SB        t
SB 1.710769 1.845449
t  1.845449 2.839359

Covariância entre os grupos: \(\mathbf{s}_{12}\)

s12<-s[1:2,3:4];s12

          SB         t
Ca 0.8628205 1.0085897
Mg 0.8261538 0.8227564

Covariância entre os grupos: \(\mathbf{s}_{21}\)

s21<-s[3:4,1:2];s21

          Ca        Mg
SB 0.8628205 0.8261538
t  1.0085897 0.8227564

Inversa da raiz quadrada de s11: \(\mathbf{s}_{11}^{-1/2}\)

library(MultBiplotR)
inv_sqrt_s11<-matrixsqrtinv(s11);inv_sqrt_s11

           [,1]       [,2]
[1,]  1.4968789 -0.4353761
[2,] -0.4353761  1.5650096

#outra alternativa
library(expm)
#inv_sqrt_s11<-solve(sqrtm(s11))

Inversa de s22: \(\mathbf{s}_{22}^{-1}\)

invs22<-solve(s22);invs22

          SB         t
SB  1.955741 -1.271139
t  -1.271139  1.178372

Matriz: \(\mathbf{s}_{11}^{-1/2}\mathbf{s}_{12}\mathbf{s}_{22}^{-1}\mathbf{s}_{21}\mathbf{s}_{11}^{-1/2}\)

matriz<-inv_sqrt_s11%*%s12%*%invs22%*%s21%*%inv_sqrt_s11;matriz

          [,1]      [,2]
[1,] 0.5328036 0.4753333
[2,] 0.4753333 0.5152602

Autovalores

Lambda_hat<-diag(eigen(matriz)$values);Lambda_hat

          [,1]       [,2]
[1,] 0.9994461 0.00000000
[2,] 0.0000000 0.04861766

Autovetores

P_hat<-eigen(matriz)$vectors;P_hat

           [,1]       [,2]
[1,] -0.7136002  0.7005531
[2,] -0.7005531 -0.7136002

Assim, as combinações lineares que definirão as variáveis canônicas do primeiro grupo de variáveis

\(\hat{C}^\top = \hat{P}\)

C_hat_t = P_hat

\(\hat{A}^\top = \mathbf{s}_{11}^{-1/2} \hat{C}^\top\)

A_hat_t <- inv_sqrt_s11 %*% C_hat_t;A_hat_t

           [,1]      [,2]
[1,] -0.7631691  1.359328
[2,] -0.7856878 -1.421795

Matriz diagonal de correlação amostral entre os dois pares de variáveis canônicas

r_U_V <- sqrt(Lambda_hat);r_U_V

         [,1]      [,2]
[1,] 0.999723 0.0000000
[2,] 0.000000 0.2204941

As correlações amostrais entre os pares \((\hat{U}_1,\hat{V}_1)\) e \((\hat{U}_2,\hat{V}_2)\) são dadas por \[r_{\hat{U}_1,\hat{V}_1} = 0.999723 \] e \[ r_{\hat{U}_2,\hat{V}_2} = 0.2204941\]

Combinações lineares para o segundo grupo

\(\hat{B}^\top = \mathbf{S}_{22}^{-1/2} \mathbf{s}_{21} \hat{A}^\top r^{-1}_{\hat{U},\hat{V}}\)

B_hat_t<- invs22%*%s21%*%A_hat_t %*% solve(r_U_V);B_hat_t 

           [,1]      [,2]
SB -0.757363058 -1.175646
t  -0.006646535  1.085508

Primeiro par de variáveis canônicas (U1,V1)

sendo \(r_{U_1,V_1} = 0.999723\)

Segundo par de variáveis canônicas

sendo \(r_{U_2,V_2} = 0.2204941\)

Correlação entre variáveis originais do primeiro grupo e variáveis canônicas de cada grupo

D11<-diag(diag(s11))
D22<-diag(diag(s22))

r_U_Y1<- t(A_hat_t)%*%s11%*%matrixsqrtinv(D11);r_U_Y1

           [,1]       [,2]
[1,] -0.8752523 -0.8719735
[2,]  0.4836667 -0.4895531

r_U_Y2<- t(A_hat_t)%*%s12%*%matrixsqrtinv(D22);r_U_Y2

             [,1]       [,2]
[1,] -0.999704256 -0.8404278
[2,] -0.001350054  0.1194113

r_V_Y2<- t(B_hat_t)%*%s22%*%matrixsqrtinv(D22);r_V_Y2

             [,1]       [,2]
[1,] -0.999981255 -0.8406606
[2,] -0.006122856  0.5415623

r_V_Y1<- t(B_hat_t)%*%s21%*%matrixsqrtinv(D11);r_V_Y1

           [,1]       [,2]
[1,] -0.8750098 -0.8717319
[2,]  0.1066457 -0.1079436

res<-data.frame(c("Ca","Mg","SB","t"),
           rbind(t(r_U_Y1),t(r_U_Y2)),rbind(t(r_V_Y1),t(r_V_Y2)))
colnames(res)<-c("Var.Originais","U_hat_1","U_hat_2","V_hat_1","V_hat_2")
res

  Var.Originais    U_hat_1      U_hat_2    V_hat_1      V_hat_2
1            Ca -0.8752523  0.483666705 -0.8750098  0.106645670
2            Mg -0.8719735 -0.489553146 -0.8717319 -0.107943595
3            SB -0.9997043 -0.001350054 -0.9999813 -0.006122856
4             t -0.8404278  0.119411304 -0.8406606  0.541562281

Exemplo 12.2

Referência: Daniel Ferreira Furtado 3a. Ed. Editora UFLA

Determinar a qualidade dos resultados da análise de variáveis canônicas

C_hat_t_1 <- C_hat_t[,1]

Lambda_hat_sqrt_1<-sqrt(Lambda_hat[1,1])

#matriz2<- s11^{-1/2} s12 s22^{-1/2} s21 s11^{-1/2}

matriz2<- inv_sqrt_s11%*%s12%*%solve(s22)%*%s21%*%inv_sqrt_s11;matriz2

          [,1]      [,2]
[1,] 0.5328036 0.4753333
[2,] 0.4753333 0.5152602

C_hat<-eigen(matriz2)$vectors

A_hat<- t(C_hat)%*%inv_sqrt_s11
  
inv_sqrt_s22 <-matrixsqrtinv(s22)

#matriz3<- s22^{-1/2} s21 s11^{-1} s12 s22^{-1/2}
matriz3<- inv_sqrt_s22%*%s21%*%solve(s11)%*%s12%*%inv_sqrt_s22;matriz3

          [,1]      [,2]
[1,] 0.7222068 0.4321404
[2,] 0.4321404 0.3258569

D_hat<-eigen(matriz3)$vectors

#k=1

A_hat_k_menos <-sqrtm(s11)%*%C_hat[,1]
B_hat_k_menos <-sqrtm(s22)%*%D_hat[,1]

#S11_til<-A_hat^{(k)-}(A_hat^{(k)-})^T = 
    #S11^{1/2}C_hat^{(k)}C_hat_t^{(k)}S11^{(1/2)}
s11_til<-A_hat_k_menos%*%t(A_hat_k_menos);s11_til

          [,1]      [,2]
[1,] 0.4360686 0.4169096
[2,] 0.4169096 0.3985924

#S22_til<-B_hat^{(k)-}(B_hat^{(k)-})^T
    #S22^{1/2}D_hat^{(k)}D_hat_t^{(k)}S22^{(1/2)}
s22_til<-B_hat_k_menos%*%t(B_hat_k_menos);s22_til

         [,1]     [,2]
[1,] 1.710705 1.852757
[2,] 1.852757 2.006604

#S12_til<- A_hat^{(k)-} sqrt(Lambda_hat^{(k)})(B_hat^{(k)-})^T
    #S11^{1/2}C_hat^{(k)}Lambda_hat^{(1/2)}D_hat_t^{(k)}S22^{(1/2)}
s12_til<-A_hat_k_menos%*%Lambda_hat_sqrt_1%*%t(B_hat_k_menos);s12_til

          [,1]      [,2]
[1,] 0.8634649 0.9351644
[2,] 0.8255279 0.8940773

Resíduos

E11 = s11-s11_til
E22 = s22-s22_til
E12 = s12-s12_til

Somas de quadrados

#SQ_E11 <-traço(t(E11)%*%E11)
SQ_E11 <-sum(diag(t(E11)%*%E11));SQ_E11

[1] 0.06697769

#SQ_E22 <-traço(t(E22)%*%E22)
SQ_E22 <-sum(diag(t(E22)%*%E22));SQ_E22

[1] 0.6935873

#SQ_E12 <-traço(t(E12)%*%E12)
SQ_E12 <-sum(diag(t(E12)%*%E12));SQ_E12

[1] 0.01047876

Proporção da explicação da variação do modelo reduzido (\(k=1\)) para a variação total de cada grupo, a partir das variáveis originais

\[R_{(1)k}^2 = \frac{\sum\limits_{j=1}^k\hat{\boldsymbol{a}}^{j\top}\hat{\boldsymbol{a}}^{j}}{tr(\boldsymbol{S}_{11})} = \frac{\sum\limits_{j=1}^k\hat{\boldsymbol{a}}^{j\top}\hat{\boldsymbol{a}}^{j}}{ \sum\limits_{j=1}^r\hat{\boldsymbol{a}}^{j\top}\hat{\boldsymbol{a}}^{j} }\]

\[R_{(2)k}^2 = \frac{\sum\limits_{j=1}^k\hat{\boldsymbol{b}}^{j\top}\hat{\boldsymbol{b}}^{j}}{tr(\boldsymbol{S}_{22})} = \frac{\sum\limits_{j=1}^k\hat{\boldsymbol{b}}^{j\top}\hat{\boldsymbol{b}}^{j}}{ \sum\limits_{j=1}^r\hat{\boldsymbol{b}}^{j\top}\hat{\boldsymbol{b}}^{j} }\] em que \(\hat{\mathbf{a}}^j = S_{11}^{1/2}\hat{\mathbf{c}}^{j\top}\) e \(\hat{\mathbf{b}}^j = S_{22}^{1/2}\hat{\mathbf{d}}^{j\top}\), com \(\hat{\mathbf{c}}^{j\top}\) e \(\hat{\mathbf{d}}^{j\top}\) os vetores linhas de \(\mathbf{C}^{k}\) e \(\mathbf{D}^{k}\), respectivamente.

tr_s11<-sum(diag(s11))
tr_s22<-sum(diag(s22))

#aj_t<-sqrtm(s11)%*%(A_hat_k_menos)
aj_t<-A_hat_k_menos
R2_1_k <-(t(aj_t)%*%aj_t)/sum(diag(s11));R2_1_k

        [,1]
[1,] 0.76332

bj_t<-B_hat_k_menos
R2_2_k<-(t(bj_t)%*%bj_t)/sum(diag(s22));R2_2_k

         [,1]
[1,] 0.816968

Proporção da explicação da variação do modelo reduzido (\(k=1\)) para a variação total de cada grupo, a partir das variáveis originais padronizadas

\[R_{(1)k}^2 = \frac{\sum\limits_{j=1}^k\sum\limits_{i=1}^p r_{\hat{U}_j, Y_{(1)i}}^2 }{ p } \]

\[R_{(2)k}^2 = \frac{\sum\limits_{j=1}^k\sum\limits_{i=1}^q r_{\hat{V}_j, Y_{(2)i}}^2 }{ q } \] \(R_{(1)k}^2\) e \(R_{(2)k}^2\) medem a importância relativa do modelo reduzido de variáveis canônicas na reprodução da covariância amostral dos grupos \(1\) e \(2\), respectivamente.

tr_s11<-sum(diag(s11))
tr_s22<-sum(diag(s22))

#aj_t<-sqrtm(s11)%*%(A_hat_k_menos)
aj_t<-A_hat_k_menos
R2_1_k <-(t(aj_t)%*%aj_t)/sum(diag(s11));R2_1_k

        [,1]
[1,] 0.76332

bj_t<-B_hat_k_menos
R2_2_k<-(t(bj_t)%*%bj_t)/sum(diag(s22));R2_2_k

         [,1]
[1,] 0.816968

Exemplo 12.3

Referência: Daniel Ferreira Furtado 3a. Ed. Editora UFLA

Aplicar o teste de nulidade dos \(r=2\) e o único teste sequencial possível se a hipótese nula for rejeitada

Inferências

Assumindo normalidade multivariada dos dados, temos

r_U_V <- sqrt(Lambda_hat);r_U_V

         [,1]      [,2]
[1,] 0.999723 0.0000000
[2,] 0.000000 0.2204941

A hipótese de interesse é dada por \[H_0:\boldsymbol{\Sigma}_{12} = \boldsymbol{0},\] que pode ser reescrita em termos das correlações canônicas, \[H_0: \rho_{U_1,V_1} = \rho_{U_2,V_2} = 0\] A estatística de teste com a correção de Bartlett é dada por \[\chi_c^2=-[n-1-(p+q+1)/2]\sum_{i=1}^{r}\ln(1-r_{\hat{U}_i,\hat{V}_i}^2)\] sendo que \(\chi_c^2\) tem distribuição aproximadamente qui-quadrado com \(\nu = pq\) graus de liberdade

chi2c<- -(n-1-(p+q+1)/2)*(log(1-r_U_V[1,1]^2)+log(1-r_U_V[2,2]^2));chi2c 

[1] 71.70894

valor_p1<-pchisq(chi2c,df=p*q,lower.tail = FALSE);valor_p1

[1] 9.887633e-15

Teste sequencial

\[H_0^1: \rho_{U_1,V_1} \neq 0, \rho_{U_2,V_2} = 0\]

\[\chi_k^2=-[n-1-(p+q+1)/2]\sum_{i=k+1}^{r}\ln(1-r_{\hat{U}_i,\hat{V}_i}^2)\]

que possui distribuição assintótica qui-quadrado com \(\nu = (p-1)(q-1)\) graus de liberdade

chi2k<- -(n-1-(p+q+1)/2)*(log(1-r_U_V[2,2]^2));chi2k 

[1] 0.473473

valor_p2<-pchisq(chi2k,df=(p-1)*(q-1),lower.tail = FALSE);valor_p2

[1] 0.4913938