Análise de Correlação Canônica

Case 1:

Dados referentes a três variáveis psicológicas (locus_of_control, self_concept, motivation), quatro variáveis acadêmicas (read, write, math, science) e uma varáveis demográfica (female) foram coletados de 600 calouros da UNICAMP. O SAE (serviço de apoio ao estudante) está interessado em saber como as variáveis psicológicas se relacionam com as variáveis acadêmicas.

O data set está disponível aqui mas está em formado *.dta (STATA) 😞.

library(dplyr)
library(haven) 
library(corrplot)
library(CCA)

# Importar data set
unicamp <- read_dta("datasets/mmreg.dta")
glimpse(unicamp)

Rows: 600
Columns: 9
$ id               <dbl> 303, 404, 225, 553, 433, 189, 336, 321, 358,…
$ locus_of_control <dbl> -0.84, -0.38, 0.89, 0.71, -0.64, 1.11, 0.06,…
$ self_concept     <dbl> -0.24, -0.47, 0.59, 0.28, 0.03, 0.90, 0.03, …
$ motivation       <dbl> 1.00, 0.67, 0.67, 0.67, 1.00, 0.33, 0.67, 0.…
$ read             <dbl> 54.8, 62.7, 60.6, 62.7, 41.6, 62.7, 41.6, 44…
$ write            <dbl> 64.5, 43.7, 56.7, 56.7, 46.3, 64.5, 39.1, 39…
$ math             <dbl> 44.5, 44.7, 70.5, 54.7, 38.4, 61.4, 56.3, 46…
$ science          <dbl> 52.6, 52.6, 58.0, 58.0, 36.3, 58.0, 45.0, 36…
$ female           <dbl> 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1,…

# Separamos os conjuntos de dados em X e Y
X <- unicamp %>% dplyr::select(locus_of_control, self_concept, motivation)
Y <- unicamp %>% dplyr::select(read, write, math, science)

# Como passo prévio, veremos se existe correlação entre as variáveis de ambos os conjuntos
unicamp %>% 
  dplyr::select(-id, -female) %>% 
  cor() %>% 
  corrplot(method = "color")

# Aplicamos CC
cc_resultados <- cc(X, Y)

# Correlação das variáveis canônicas
cc_resultados$cor

[1] 0.44643648 0.15335902 0.02250348

• As primeiras variávais canônicas (\(\eta_1\) e \(\phi_1\)) tem uma correlação de 0.45.
• As segundas variáveis canônicas (\(\eta_2\) e \(\phi_2\)) tem uma correlação de 0.15.
• As terceiras variáveis canônicas (\(\eta_3\) e \(\phi_3\)) tem uma correlação de 0.02.

# Quais as combinações lineares para formar as variáves canónicas?
cc_resultados$xcoef

                       [,1]       [,2]       [,3]
locus_of_control -1.2501212  0.7659633  0.4966529
self_concept      0.2367331  0.8421110 -1.2051226
motivation       -1.2491434 -2.6359626 -1.0935084

• A uma unidade de aumento em locus_of_control a primeira variável canônica diminui em 1.25 pontos.
• A uma unidade de aumento em motivation a primeira variável canônica diminui em 1.249 pontos.
• A uma unidade de aumento em self_concept a primeira variável canônica aumenta em 0.24 pontos.
• A primeira variável canônica é formada, principalmente por, locus_of_control e motivation.

# Quais as combinações lineares para formar as variáves canónicas?
cc_resultados$ycoef

                [,1]         [,2]        [,3]
read    -0.044047126 -0.001592913 -0.08833162
write   -0.055088843 -0.090414603  0.09612884
math    -0.019401099 -0.002955461 -0.08782254
science  0.003797757  0.124208982  0.08849518

• Quando read aumenta em uma unidade, a primeira variável canônica diminui em 0.044 unidades.
• Quando write aumenta em uma unidade, a primeira variável canônica diminui em 0.055 unidades.
• Quando math aumenta em uma unidade, a primeira variável canônica diminui em 0.019 unidades.

Até aqui, temos feito tudo o visto nos slides. Contudo, sob normalidade multivariada, podemos ir além e fazer inferência estatística.

Testar se X e Y são não correlacionados.

• Sob normalidade e utilizando o TRV chegamos na seguinte estatística de teste \[|\textbf{I} - \textbf{S}_{YY}^{-1}\textbf{S}_{YX} \textbf{S}_{XX}^{-1} \textbf{SY}| = \displaystyle \prod_{i = 1}^k(1-l_i).\]
• Utilizando a aproximação de Barlett (1939), temos que \[-[n - (p + q + 3)/2] \displaystyle \sum_{i = 1}^k \log(1-l_i) \sim \chi^2_{pq}\]

mvnormalTest::mardia(cbind(X,Y))

$mv.test

$uv.shapiro
                 W      p-value UV.Normality
locus_of_control 0.971  0       No          
self_concept     0.941  0       No          
motivation       0.8238 0       No          
read             0.9816 0       No          
write            0.9533 0       No          
math             0.9835 0       No          
science          0.983  0       No          

aux <- as.matrix(cbind(X,Y))
mnt::test.BHEP(aux, MC.rep = 1000)  # Rejeitamos H0

Os dados não são normalmente distribuidos. Apenas a maneira de exemplo continuaremos com inferência estatística.

n <- nrow(unicamp)
p <- ncol(X)
q <- ncol(Y)
l <- cc_resultados$cor^2 #(Lembre que SVD fornece em D a raiz dos autovalores)
e_teste <- -(n - (p + q + 3)/2)*sum(log(1 - l))
e_teste

[1] 146.7165

x2 <- qchisq(0.95, p*q)
ifelse(e_teste > x2, "Rejeitamos H0", "Não rejeitamos H0")

[1] "Rejeitamos H0"

Testar que a correlaão de apenas k - s variáveis canônicas são diferentes de zero.

\[-[n - (p + q + 3)/2] \displaystyle \sum_{i = s + 1}^k \log(1-l_i) \sim \chi^2_{(p - s)(q - s)}\]

s <- 0  # Caso anterior
k <- 3
e_teste <- -(n - (p + q + 3)/2)*sum(log(1 - l)[(s + 1):k])
e_teste

[1] 146.7165

x2 <- qchisq(0.95, (p - s)*(q - s))
ifelse(e_teste > x2, "Rejeitamos H0", "Não rejeitamos H0")

[1] "Rejeitamos H0"

s <- 1  
k <- 3
e_teste <- -(n - (p + q + 3)/2)*sum(log(1 - l)[(s + 1):k])
e_teste

[1] 14.46237

x2 <- qchisq(0.95, (p - s)*(q - s))
ifelse(e_teste > x2, "Rejeitamos H0", "Não rejeitamos H0")

[1] "Rejeitamos H0"

s <- 2  
k <- 3
e_teste <- -(n - (p + q + 3)/2)*sum(log(1 - l)[(s + 1):k])
e_teste

[1] 0.3013883

x2 <- qchisq(0.95, (p - s)*(q - s))
ifelse(e_teste > x2, "Rejeitamos H0", "Não rejeitamos H0")

[1] "Não rejeitamos H0"

• Apenas as 2 primeiras variáveis canônicas tem correlação diferentes de zero.

Case 2:

Roberto Fulcherberguer, o CEO das Casas Bahia, está interessado em saber se existe alguma relação entre o desempenho das vendas (mensurado através das variáveis Sales Growth, Sales Profitability, New Account Sales) e a inteligência dos gerentes da loja (mensurado através das variáveis Creativity, Mechanical Reasoning, Abstract Reasoning, Mathematics). Para fazer as análises, \(n = 50\) lojas são selecionadas aleatoriamente e as variáveis de interesse são obtidas.

O data set está disponível aqui.

# Importar dados
casas_bahia <- read.table("datasets/sales.txt")
colnames(casas_bahia) <- c("growth", "profit", "new",
                           "create", "mechanical", "abs", "math")
glimpse(casas_bahia)

Rows: 50
Columns: 7
$ growth     <dbl> 93.0, 88.8, 95.0, 101.3, 1…
$ profit     <dbl> 96.0, 91.8, 100.3, 103.8, …
$ new        <dbl> 97.8, 96.8, 99.0, 106.8, 1…
$ create     <int> 9, 7, 8, 13, 10, 10, 9, 18…
$ mechanical <int> 12, 10, 12, 14, 15, 14, 12…
$ abs        <int> 9, 10, 9, 12, 12, 11, 9, 1…
$ math       <int> 20, 15, 26, 29, 32, 21, 25…

# Separamos os conjuntos de dados em X e Y
X <- casas_bahia %>% dplyr::select(growth, profit, new)
Y <- casas_bahia %>% dplyr::select(create, mechanical, abs, math)

# Correlação entre as variáveis
casas_bahia %>% 
  cor() %>% 
  corrplot(method = "color")

# Aplicamos CC
cc_cb <- cc(X, Y)

# Correlação das variáveis canônicas
cc_cb$cor

[1] 0.9944827 0.8781065 0.3836057

• As primeiras variáveis canônicas tem uma correlação de 0.99.
• As segundas variáveis canônicas tem uma correlação de 0.88.
• As terceiras variáveis canônicas tem uma correlação de 0.28.

# Quais as combinações lineares para formar as variáves canónicas?
cc_cb$xcoef

              [,1]       [,2]       [,3]
growth -0.06237788 -0.1740703  0.3771529
profit -0.02092564  0.2421641 -0.1035150
new    -0.07825817 -0.2382940 -0.3834151

# Quais as combinações lineares para formar as variáves canónicas?
cc_cb$ycoef

                  [,1]        [,2]        [,3]
create     -0.06974814 -0.19239132 -0.24655659
mechanical -0.03073830  0.20157438  0.14189528
abs        -0.08956418 -0.49576326  0.28022405
math       -0.06282997  0.06831607 -0.01133259

mvnormalTest::mardia(cbind(X,Y))

$mv.test

$uv.shapiro
           W      p-value UV.Normality
growth     0.9401 0.0135  No          
profit     0.9557 0.0589  Yes         
new        0.9781 0.4732  Yes         
create     0.9604 0.0921  Yes         
mechanical 0.9723 0.287   Yes         
abs        0.9482 0.0287  No          
math       0.9768 0.4255  Yes