Tema: Correlação
Prof Dr. Roberto C. Leoni
14 de outubro de 2016
Este material apresenta vários exemplos comentados para calcular correlação de Pearson, Spearman, Kendal, teste de correlação, teste de normalidade e gráficos diversos sobre o tema.
Não é meu objetivo apresentar a parte teórica sobre o tema. Busque outras fontes para melhor compeensão das técnicas apresentadas a seguir.
instale os pacotes a seguir a primeira vex que for executar os exemplos
# install.packages("ggpubr","Hmisc", "corrplot","PerformanceAnalytics") # carregar o pacote
library("ggpubr")Vamos usar o pacote mtcars para ser a nossa base de dados. São 11 caractéristicas de 32 veículos
# .txt tab file, use read.delim
# my_data <- read.delim(file.choose())
# .csv file, use read.csv
# my_data <- read.csv(file.choose())
my_data <- mtcars
head(my_data, 6)## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1O gráfico a seguir representa a relação entre as variáveis mpg e wt (apresenta a correlação de pearson e o p-valor). Explicando de forma prática o p-valor: Quando é próximo a zero, significa que não há evidência estatística de que a correlação entre as variáveis seja igual a zero. Quanto maior for o p-valor, menor é a chance de evidenciar essa hipótese e vice-versa.
ggscatter(my_data, x = "mpg", y = "wt",
add = "reg.line", conf.int = TRUE,
cor.coef = TRUE, cor.method = "pearson",
xlab = "Miles/(US) gallon", ylab = "Weight (1000 lbs)")
o coeficiente de correlação de Pearson exige que os dados sigam a distribuição normal de probabilidades. O teste a seguir (Shapiro Wilk) é uma técnica formal para avaliar normalidade. Se o p-valor for próximo a zero, rejeitamos a hipótese de normalidade. Um valor < 0,05, por exemplo, geralmente é empregado como referência.
# Shapiro-Wilk normality test for mpg
shapiro.test(my_data$mpg) # => p = 0.1229##
## Shapiro-Wilk normality test
##
## data: my_data$mpg
## W = 0.94756, p-value = 0.1229# Shapiro-Wilk normality test for wt
shapiro.test(my_data$wt) # => p = 0.09##
## Shapiro-Wilk normality test
##
## data: my_data$wt
## W = 0.94326, p-value = 0.09265Uma forma empírica para avaliar normalidade é o gráfico qqplot. Haverá evidência empírica de normalidade na medida em que os valores estejam próximmos da reta.
# mpg
ggqqplot(my_data, my_data$mpg, ylab = "MPG")
# wt
ggqqplot(my_data, my_data$wt, ylab = "WT")
teste da correlação de Pearson. Apresenta a correlação e o p-valor
res <- cor.test(my_data$wt, my_data$mpg,
method = "pearson")
res##
## Pearson's product-moment correlation
##
## data: my_data$wt and my_data$mpg
## t = -9.559, df = 30, p-value = 1.294e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9338264 -0.7440872
## sample estimates:
## cor
## -0.8676594Os testes de kendal e Spearman são medidas de correlações emregadas para dados ordinais (Ex: escalas de likert de 1 a 5, ordens de classificação, julgamentos realizados em categorias e etc. A ideia é que os resultados representem hierarquia, ou seja, 1 representa um resultado inferior a 2 e assim por diante. Contudo não podemos afirmar que o 2 é duas vezes melhor que o 1), ou ainda, estes testes são epregados quando não podemos supor normalidade nos dados.
teste de correlação de kendall. Apresenta a correlação e o p-valor
res2 <- cor.test(my_data$wt, my_data$mpg, method="kendall")
res2##
## Kendall's rank correlation tau
##
## data: my_data$wt and my_data$mpg
## z = -5.7981, p-value = 6.706e-09
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
## tau
## -0.7278321teste de correlação de spearman. Apresenta a correlação e o p-valor
res3 <-cor.test(my_data$wt, my_data$mpg, method = "spearman")
res3##
## Spearman's rank correlation rho
##
## data: my_data$wt and my_data$mpg
## S = 10292, p-value = 1.488e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.886422matriz com a correlação entre todas as variáveis
res <- cor(my_data, use = "complete.obs")
round(res, 2)## mpg cyl disp hp drat wt qsec vs am gear carb
## mpg 1.00 -0.85 -0.85 -0.78 0.68 -0.87 0.42 0.66 0.60 0.48 -0.55
## cyl -0.85 1.00 0.90 0.83 -0.70 0.78 -0.59 -0.81 -0.52 -0.49 0.53
## disp -0.85 0.90 1.00 0.79 -0.71 0.89 -0.43 -0.71 -0.59 -0.56 0.39
## hp -0.78 0.83 0.79 1.00 -0.45 0.66 -0.71 -0.72 -0.24 -0.13 0.75
## drat 0.68 -0.70 -0.71 -0.45 1.00 -0.71 0.09 0.44 0.71 0.70 -0.09
## wt -0.87 0.78 0.89 0.66 -0.71 1.00 -0.17 -0.55 -0.69 -0.58 0.43
## qsec 0.42 -0.59 -0.43 -0.71 0.09 -0.17 1.00 0.74 -0.23 -0.21 -0.66
## vs 0.66 -0.81 -0.71 -0.72 0.44 -0.55 0.74 1.00 0.17 0.21 -0.57
## am 0.60 -0.52 -0.59 -0.24 0.71 -0.69 -0.23 0.17 1.00 0.79 0.06
## gear 0.48 -0.49 -0.56 -0.13 0.70 -0.58 -0.21 0.21 0.79 1.00 0.27
## carb -0.55 0.53 0.39 0.75 -0.09 0.43 -0.66 -0.57 0.06 0.27 1.00matriz com a correlação entre todas as variáveis
library("Hmisc")
res2 <- rcorr(as.matrix(my_data))
res2## mpg cyl disp hp drat wt qsec vs am gear carb
## mpg 1.00 -0.85 -0.85 -0.78 0.68 -0.87 0.42 0.66 0.60 0.48 -0.55
## cyl -0.85 1.00 0.90 0.83 -0.70 0.78 -0.59 -0.81 -0.52 -0.49 0.53
## disp -0.85 0.90 1.00 0.79 -0.71 0.89 -0.43 -0.71 -0.59 -0.56 0.39
## hp -0.78 0.83 0.79 1.00 -0.45 0.66 -0.71 -0.72 -0.24 -0.13 0.75
## drat 0.68 -0.70 -0.71 -0.45 1.00 -0.71 0.09 0.44 0.71 0.70 -0.09
## wt -0.87 0.78 0.89 0.66 -0.71 1.00 -0.17 -0.55 -0.69 -0.58 0.43
## qsec 0.42 -0.59 -0.43 -0.71 0.09 -0.17 1.00 0.74 -0.23 -0.21 -0.66
## vs 0.66 -0.81 -0.71 -0.72 0.44 -0.55 0.74 1.00 0.17 0.21 -0.57
## am 0.60 -0.52 -0.59 -0.24 0.71 -0.69 -0.23 0.17 1.00 0.79 0.06
## gear 0.48 -0.49 -0.56 -0.13 0.70 -0.58 -0.21 0.21 0.79 1.00 0.27
## carb -0.55 0.53 0.39 0.75 -0.09 0.43 -0.66 -0.57 0.06 0.27 1.00
##
## n= 32
##
##
## P
## mpg cyl disp hp drat wt qsec vs am gear
## mpg 0.0000 0.0000 0.0000 0.0000 0.0000 0.0171 0.0000 0.0003 0.0054
## cyl 0.0000 0.0000 0.0000 0.0000 0.0000 0.0004 0.0000 0.0022 0.0042
## disp 0.0000 0.0000 0.0000 0.0000 0.0000 0.0131 0.0000 0.0004 0.0010
## hp 0.0000 0.0000 0.0000 0.0100 0.0000 0.0000 0.0000 0.1798 0.4930
## drat 0.0000 0.0000 0.0000 0.0100 0.0000 0.6196 0.0117 0.0000 0.0000
## wt 0.0000 0.0000 0.0000 0.0000 0.0000 0.3389 0.0010 0.0000 0.0005
## qsec 0.0171 0.0004 0.0131 0.0000 0.6196 0.3389 0.0000 0.2057 0.2425
## vs 0.0000 0.0000 0.0000 0.0000 0.0117 0.0010 0.0000 0.3570 0.2579
## am 0.0003 0.0022 0.0004 0.1798 0.0000 0.0000 0.2057 0.3570 0.0000
## gear 0.0054 0.0042 0.0010 0.4930 0.0000 0.0005 0.2425 0.2579 0.0000
## carb 0.0011 0.0019 0.0253 0.0000 0.6212 0.0146 0.0000 0.0007 0.7545 0.1290
## carb
## mpg 0.0011
## cyl 0.0019
## disp 0.0253
## hp 0.0000
## drat 0.6212
## wt 0.0146
## qsec 0.0000
## vs 0.0007
## am 0.7545
## gear 0.1290
## carbp-valores e Correlações
# Extract the correlation coefficients
res2$r## mpg cyl disp hp drat wt
## mpg 1.0000000 -0.8521619 -0.8475513 -0.7761683 0.68117189 -0.8676594
## cyl -0.8521619 1.0000000 0.9020329 0.8324475 -0.69993812 0.7824958
## disp -0.8475513 0.9020329 1.0000000 0.7909486 -0.71021390 0.8879799
## hp -0.7761683 0.8324475 0.7909486 1.0000000 -0.44875914 0.6587479
## drat 0.6811719 -0.6999381 -0.7102139 -0.4487591 1.00000000 -0.7124406
## wt -0.8676594 0.7824958 0.8879799 0.6587479 -0.71244061 1.0000000
## qsec 0.4186840 -0.5912421 -0.4336979 -0.7082234 0.09120482 -0.1747159
## vs 0.6640389 -0.8108118 -0.7104159 -0.7230967 0.44027850 -0.5549157
## am 0.5998324 -0.5226070 -0.5912271 -0.2432043 0.71271110 -0.6924953
## gear 0.4802848 -0.4926866 -0.5555692 -0.1257043 0.69961011 -0.5832870
## carb -0.5509251 0.5269883 0.3949769 0.7498125 -0.09078978 0.4276059
## qsec vs am gear carb
## mpg 0.41868404 0.6640389 0.59983242 0.4802848 -0.55092508
## cyl -0.59124213 -0.8108118 -0.52260703 -0.4926866 0.52698827
## disp -0.43369791 -0.7104159 -0.59122705 -0.5555692 0.39497685
## hp -0.70822340 -0.7230967 -0.24320425 -0.1257043 0.74981248
## drat 0.09120482 0.4402785 0.71271110 0.6996101 -0.09078978
## wt -0.17471591 -0.5549157 -0.69249529 -0.5832870 0.42760593
## qsec 1.00000000 0.7445354 -0.22986083 -0.2126822 -0.65624928
## vs 0.74453545 1.0000000 0.16834512 0.2060234 -0.56960714
## am -0.22986083 0.1683451 1.00000000 0.7940587 0.05753435
## gear -0.21268222 0.2060234 0.79405874 1.0000000 0.27407283
## carb -0.65624928 -0.5696071 0.05753435 0.2740728 1.00000000# Extract p-values
res2$P## mpg cyl disp hp drat
## mpg NA 6.112697e-10 9.380354e-10 1.787838e-07 1.776241e-05
## cyl 6.112697e-10 NA 1.803002e-12 3.477856e-09 8.244635e-06
## disp 9.380354e-10 1.803002e-12 NA 7.142686e-08 5.282028e-06
## hp 1.787838e-07 3.477856e-09 7.142686e-08 NA 9.988768e-03
## drat 1.776241e-05 8.244635e-06 5.282028e-06 9.988768e-03 NA
## wt 1.293956e-10 1.217567e-07 1.222311e-11 4.145833e-05 4.784268e-06
## qsec 1.708199e-02 3.660527e-04 1.314403e-02 5.766250e-06 6.195823e-01
## vs 3.415940e-05 1.843015e-08 5.235010e-06 2.940897e-06 1.167552e-02
## am 2.850209e-04 2.151208e-03 3.662112e-04 1.798309e-01 4.726797e-06
## gear 5.400949e-03 4.173297e-03 9.635927e-04 4.930119e-01 8.360117e-06
## carb 1.084446e-03 1.942341e-03 2.526789e-02 7.827805e-07 6.211834e-01
## wt qsec vs am gear
## mpg 1.293956e-10 1.708199e-02 3.415940e-05 2.850209e-04 5.400949e-03
## cyl 1.217567e-07 3.660527e-04 1.843015e-08 2.151208e-03 4.173297e-03
## disp 1.222311e-11 1.314403e-02 5.235010e-06 3.662112e-04 9.635927e-04
## hp 4.145833e-05 5.766250e-06 2.940897e-06 1.798309e-01 4.930119e-01
## drat 4.784268e-06 6.195823e-01 1.167552e-02 4.726797e-06 8.360117e-06
## wt NA 3.388682e-01 9.798495e-04 1.125438e-05 4.586600e-04
## qsec 3.388682e-01 NA 1.029669e-06 2.056622e-01 2.425345e-01
## vs 9.798495e-04 1.029669e-06 NA 3.570440e-01 2.579439e-01
## am 1.125438e-05 2.056622e-01 3.570440e-01 NA 5.834051e-08
## gear 4.586600e-04 2.425345e-01 2.579439e-01 5.834051e-08 NA
## carb 1.463861e-02 4.536940e-05 6.670496e-04 7.544526e-01 1.290291e-01
## carb
## mpg 1.084446e-03
## cyl 1.942341e-03
## disp 2.526789e-02
## hp 7.827805e-07
## drat 6.211834e-01
## wt 1.463861e-02
## qsec 4.536940e-05
## vs 6.670496e-04
## am 7.544526e-01
## gear 1.290291e-01
## carb NAa função a seguir flattenCorrMatrix será empregada para organizar a saída.
# ++++++++++++++++++++++++++++
# flattenCorrMatrix
# ++++++++++++++++++++++++++++
# cormat : matrix of the correlation coefficients
# pmat : matrix of the correlation p-values
flattenCorrMatrix <- function(cormat, pmat) {
ut <- upper.tri(cormat)
data.frame(
row = rownames(cormat)[row(cormat)[ut]],
column = rownames(cormat)[col(cormat)[ut]],
cor =(cormat)[ut],
p = pmat[ut]
)
}Emprego da função ** flattenCorrMatrix **
library(Hmisc)
res2<-rcorr(as.matrix(mtcars[,1:7])) # somente as variáveis 1 a 7
flattenCorrMatrix(res2$r, res2$P)## row column cor p
## 1 mpg cyl -0.85216194 6.112697e-10
## 2 mpg disp -0.84755135 9.380354e-10
## 3 cyl disp 0.90203285 1.803002e-12
## 4 mpg hp -0.77616835 1.787838e-07
## 5 cyl hp 0.83244747 3.477856e-09
## 6 disp hp 0.79094857 7.142686e-08
## 7 mpg drat 0.68117189 1.776241e-05
## 8 cyl drat -0.69993812 8.244635e-06
## 9 disp drat -0.71021390 5.282028e-06
## 10 hp drat -0.44875914 9.988768e-03
## 11 mpg wt -0.86765939 1.293956e-10
## 12 cyl wt 0.78249580 1.217567e-07
## 13 disp wt 0.88797992 1.222311e-11
## 14 hp wt 0.65874785 4.145833e-05
## 15 drat wt -0.71244061 4.784268e-06
## 16 mpg qsec 0.41868404 1.708199e-02
## 17 cyl qsec -0.59124213 3.660527e-04
## 18 disp qsec -0.43369791 1.314403e-02
## 19 hp qsec -0.70822340 5.766250e-06
## 20 drat qsec 0.09120482 6.195823e-01
## 21 wt qsec -0.17471591 3.388682e-01Gráfico que ilustra as correlações corplot
library(corrplot)
corrplot(res, type = "upper", order = "hclust",
tl.col = "black", tl.srt = 45)
Gráfico que ilustra as correlações corplot. Mantém apenas as correlções significativas. o ponto de corte é p-valor maior que 0,01.
# Insignificant correlation are crossed
corrplot(res2$r, type="upper", order="hclust",
p.mat = res2$P, sig.level = 0.01, insig = "blank")
Outra forma gráfica para ilustras as correlações
library("PerformanceAnalytics")
my_data <- mtcars[, c(1,3,4,5,6,7)]
chart.Correlation(my_data, histogram=TRUE, pch=19)
o heatmap também ilustra as correlações. Quanto mais azul, maior é a associação positiva e vermelho representa a associação negativa (inversa). A cor branca indica ausência de correlação.
# Get some colors
col<- colorRampPalette(c("blue", "white", "red"))(20)
heatmap(x = res, col = col, symm = TRUE)
Várias formas para ilustrar a correlação. Explore!
head(mtcars)## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1M <-cor(mtcars)
head(round(M,2))## mpg cyl disp hp drat wt qsec vs am gear carb
## mpg 1.00 -0.85 -0.85 -0.78 0.68 -0.87 0.42 0.66 0.60 0.48 -0.55
## cyl -0.85 1.00 0.90 0.83 -0.70 0.78 -0.59 -0.81 -0.52 -0.49 0.53
## disp -0.85 0.90 1.00 0.79 -0.71 0.89 -0.43 -0.71 -0.59 -0.56 0.39
## hp -0.78 0.83 0.79 1.00 -0.45 0.66 -0.71 -0.72 -0.24 -0.13 0.75
## drat 0.68 -0.70 -0.71 -0.45 1.00 -0.71 0.09 0.44 0.71 0.70 -0.09
## wt -0.87 0.78 0.89 0.66 -0.71 1.00 -0.17 -0.55 -0.69 -0.58 0.43library(corrplot)
corrplot(M, method="circle")
corrplot(M, method="pie")
corrplot(M, method="color")
corrplot(M, method="number")
corrplot(M, type="upper")
corrplot(M, type="lower")