Correlação

Serve para estudar o comportamento conjunto de duas variáveis quantitativas distintas. Em outras palavras, mede o grau de associação entre duas variáveis X e Y.

Gráfico de dispersão

Em primeiro lugar devemos criar um gráfico de dispersão para analisar o conjunto de dados.

Ele permite verificar se há algum comportamento conhecido entre as variáveis.

Cada ponto é chamado de coordenada, pois representa ao mesmo tempo a variável X e Y. Por exemplo o ponto A(4,4), indica que no eixo x o valor da variável é 4 e no eixo y também vale 4. O ponto c(-3,-3), indica que a variável x é -3 e a variável y é -3.

Se representa todo o conjunto de dados em pares e se faz uma análise para ver o seu comportamento.

Exemplo no R

idade=c(12,16,21,23,23,28,29)
peso=c(56,78,68,62,91,105,71)
tabela=cbind(idade,peso)
tabela
##      idade peso
## [1,]    12   56
## [2,]    16   78
## [3,]    21   68
## [4,]    23   62
## [5,]    23   91
## [6,]    28  105
## [7,]    29   71

O que devemos fazer é plotar cada par ordenado no gráfico

plot(idade,peso,col='#FF2400',lwd=7,xlim=c(8,32),ylim=c(40,120))
text(idade,peso,cex=1.6)

Há três situações que são muito marcantes de acontecer:

Se os pontos estiverem dispersos, sem definição de direção, então podemos dizer que há uma correlação fraca ou nula.

Para se medir essa associação, se calcula o Coefifiente de Pearson que pode variar de -1 ate 1. Quanto mais próximo de 1, será forte positiva; e quanto mais próximo de -1 será forte negativa.

Para se calcular usamos a seguinte fórmula:

\[r=\frac{ \sum_{i=1}^n X_i Y_i - \frac{\sum_{i=1}^n X_i \sum_{i=1}^n Y_i}{n}}{\sqrt{(\sum_{i=1}^n X_i^2 -\frac{(\sum_{i=1} X_i)^2}{n}) (\sum_{i=1}^n Y_i^2 -\frac{(\sum_{i=1} Y_i)^2}{n}) }}\]

com \(-1 \leq r \leq 1\).

Exemplo no R

Vamos fazer a correlação entre as variáveis idade e peso para verificar a correlação.

cor(idade,peso)
## [1] 0.5361417

Notamos que a correlação é positiva, porém moderada.

Mesmo assim ainda é necessário ver se há significância na correlação por teste de hipóteses.

Assim usamos o teste de correlação.

cor.test(idade,peso)
## 
##  Pearson's product-moment correlation
## 
## data:  idade and peso
## t = 1.4202, df = 5, p-value = 0.2148
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3637986  0.9183997
## sample estimates:
##       cor 
## 0.5361417

Conclusão: Como o valor-p é maior que alpha, então identificamos que não há correlação alguma, e temos \(r=0\).

if(!require(ggpubr)){install.packages("ggpubr")}
## Loading required package: ggpubr
## Loading required package: ggplot2
if(!require(Hmisc)){install.packages("Hmisc")}
## Loading required package: Hmisc
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
if(!require(corrplot)){install.packages("corrplot")}
## Loading required package: corrplot
## corrplot 0.84 loaded
if(!require(PerformanceAnalytics)){install.packages("PerformanceAnalytics")}
## Loading required package: PerformanceAnalytics
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend

Exemplo 2

Vamos avaliar o tempo(em horas) de estudo de uma turma de alunos, e suas notas obtidas em uma avaliação.

tempo=c(3,7,2,1.5,12)
nota=c(4.5,6.5,3.7,4,9.3)
boletim=cbind(tempo,nota)
boletim
##      tempo nota
## [1,]   3.0  4.5
## [2,]   7.0  6.5
## [3,]   2.0  3.7
## [4,]   1.5  4.0
## [5,]  12.0  9.3

Vamos plotar o gráfico de dispersão.

plot(tempo,nota,lwd=7,xlim=c(1,13),ylim=c(0,10),col="#FF5671",xlab = "tempo de estudo",ylab = "nota da avaliação",main=" Boletim escolar")

Notamos que quanto maior o empenho do aluno, maior a sua nota na avaliação. Assim queremos ver se há realmente uma correlação entre as duas variáveis.

cor(tempo,nota)
## [1] 0.9960249
cor.test(tempo,nota)
## 
##  Pearson's product-moment correlation
## 
## data:  tempo and nota
## t = 19.367, df = 3, p-value = 0.0003007
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9382859 0.9997509
## sample estimates:
##       cor 
## 0.9960249

Assim ela é forte e é significativa.

A análise de regressão

A análise de regressão estuda a relação entre uma variável chamada a variável dependente(Y) e outras variáveis chamadas variáveis independentes(X).

A relação entre elas é representada por um modelo matemático, que associa a variável dependente com as variáveis independentes.

Este modelo é designado por modelo de regressão linear simples se define uma relação linear entre a variável dependente e uma variável independente.

Se em vez de uma, forem incorporadas várias variáveis independentes, o modelo passa a denominar-se modelo de regressão linear múltipla.

Exemplos:

Sob dois pontos de vista:

O modelo matemático

Y - variável explicada ou dependente (aleatória);

X - variável explicativa ou independente medida sem erro (não aleatória);

\(\alpha\) ou \(b_0\) - coeficiente de regressão, que representa o intercepto (parâmetro desconhecido do modelo -> a estimar)

\(\beta\) ou \(b_1\) - coeficiente de regressão, que representa o declive (inclinação) (parâmetro desconhecido do modelo -> a estimar)

\(\epsilon\) - erro aleatório ou estocástico, onde se procuram incluir todas as influências no comportamento da variável Y que não podem ser explicadas linearmente pelo comportamento da variável X;

Os coeficientes

O cálculo dos coeficientes pode se fazer por meio:

\[\beta=\frac{\sum x_i y_i - n \bar{x} \bar{y}}{\sum x_i^2 - n {\bar{x}}^2}\]

e

\[\alpha=\bar{y}- \beta \bar{x}\]

Obtendo o modelo (modelo ajustado)

\[y_i^{prev} = \alpha + \beta x + \epsilon\]

A diferençã entre \(y_i\) e \(y_i^{prev}\) é chamada de resíduos.

\[res= y_i - y_i^{prev}\]

Também é construida uma tabela de ANOVA para o modelo de regressão, pois se precisa avaliar a significância do modelo e dos coeficientes.

Exemplo no R

modelo=lm(nota~tempo)
anova(modelo)
summary(modelo)
## 
## Call:
## lm(formula = nota ~ tempo)
## 
## Residuals:
##         1         2         3         4         5 
##  0.006394 -0.101023 -0.266752  0.296675  0.064706 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.9130     0.1756   16.59 0.000476 ***
## tempo         0.5269     0.0272   19.37 0.000301 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2406 on 3 degrees of freedom
## Multiple R-squared:  0.9921, Adjusted R-squared:  0.9894 
## F-statistic: 375.1 on 1 and 3 DF,  p-value: 0.0003007
plot(tempo,nota,col="red")
points(tempo,modelo$fitted.values,col="blue")
abline(modelo,col=5)

Avaliando os resíduos

Ao gerar os gráficos dos resíduos, podemos avaliar se o modelo é bem ajustado.

par(mfrow=c(2,2))
plot(modelo)

Exemplo do bando de dados Pizza

setwd("C:\\Users\\ednas\\OneDrive\\Documentos\\_TRABALHO\\IFMT\\MESTRADO\\Aulas_Mestrado\\Aula Red\\Aulas_latex\\regressao")

pizza=read.table("pizza.txt",h=T)

Vamos usar o banco de dados pizza que já conhecemos para avaliar uma correlação e um modelo de regressão.

head(pizza)
str(pizza)
## 'data.frame':    300 obs. of  13 variables:
##  $ marca       : chr  "A" "A" "A" "A" ...
##  $ id          : int  14069 14053 14025 14016 14005 14075 14082 14097 14117 14133 ...
##  $ agua_100    : num  27.8 28.5 28.4 30.6 30.5 ...
##  $ proteina    : num  21.4 21.3 20 20.1 21.3 ...
##  $ gordura     : num  44.9 43.9 45.8 43.1 41.6 ...
##  $ cinzas      : num  5.11 5.34 5.08 4.79 4.82 4.92 4.71 5.28 5.02 5.16 ...
##  $ sodio       : num  1.77 1.79 1.63 1.61 1.64 1.65 1.58 1.75 1.71 1.66 ...
##  $ carboidratos: num  0.77 1.02 0.8 1.38 1.76 1.4 1.77 2.95 1.18 0.64 ...
##  $ calorias    : num  4.93 4.84 4.95 4.74 4.67 4.67 4.63 4.72 4.93 4.95 ...
##  $ fds         : int  1 0 0 0 1 0 1 0 1 1 ...
##  $ fds_fator   : chr  "final_de_semana" "meio_de_semana" "meio_de_semana" "meio_de_semana" ...
##  $ borda       : int  1 1 0 0 0 1 0 0 1 1 ...
##  $ borda_fator : chr  "com_borda" "com_borda" "sem_borda" "sem_borda" ...

Precisamos montar uma tabela que tenha somente as variáveis numéricas.

tabela_pizza=pizza[,-c(1,2,11,13)]
head(tabela_pizza)
str(tabela_pizza)
## 'data.frame':    300 obs. of  9 variables:
##  $ agua_100    : num  27.8 28.5 28.4 30.6 30.5 ...
##  $ proteina    : num  21.4 21.3 20 20.1 21.3 ...
##  $ gordura     : num  44.9 43.9 45.8 43.1 41.6 ...
##  $ cinzas      : num  5.11 5.34 5.08 4.79 4.82 4.92 4.71 5.28 5.02 5.16 ...
##  $ sodio       : num  1.77 1.79 1.63 1.61 1.64 1.65 1.58 1.75 1.71 1.66 ...
##  $ carboidratos: num  0.77 1.02 0.8 1.38 1.76 1.4 1.77 2.95 1.18 0.64 ...
##  $ calorias    : num  4.93 4.84 4.95 4.74 4.67 4.67 4.63 4.72 4.93 4.95 ...
##  $ fds         : int  1 0 0 0 1 0 1 0 1 1 ...
##  $ borda       : int  1 1 0 0 0 1 0 0 1 1 ...
attach(tabela_pizza)

Agora vamos fazer uma tabela de correlação entre todas as numéricas.

round(cor(tabela_pizza),2)
##              agua_100 proteina gordura cinzas sodio carboidratos calorias   fds
## agua_100         1.00     0.36   -0.17   0.27 -0.10        -0.59    -0.76 -0.06
## proteina         0.36     1.00    0.50   0.82  0.43        -0.85     0.07 -0.05
## gordura         -0.17     0.50    1.00   0.79  0.93        -0.64     0.76 -0.04
## cinzas           0.27     0.82    0.79   1.00  0.81        -0.90     0.33 -0.09
## sodio           -0.10     0.43    0.93   0.81  1.00        -0.62     0.67 -0.07
## carboidratos    -0.59    -0.85   -0.64  -0.90 -0.62         1.00    -0.02  0.08
## calorias        -0.76     0.07    0.76   0.33  0.67        -0.02     1.00  0.02
## fds             -0.06    -0.05   -0.04  -0.09 -0.07         0.08     0.02  1.00
## borda            0.08     0.07    0.01   0.05  0.03        -0.08    -0.04  0.02
##              borda
## agua_100      0.08
## proteina      0.07
## gordura       0.01
## cinzas        0.05
## sodio         0.03
## carboidratos -0.08
## calorias     -0.04
## fds           0.02
## borda         1.00
rcorr(as.matrix(tabela_pizza))
##              agua_100 proteina gordura cinzas sodio carboidratos calorias   fds
## agua_100         1.00     0.36   -0.17   0.27 -0.10        -0.59    -0.76 -0.06
## proteina         0.36     1.00    0.50   0.82  0.43        -0.85     0.07 -0.05
## gordura         -0.17     0.50    1.00   0.79  0.93        -0.64     0.76 -0.04
## cinzas           0.27     0.82    0.79   1.00  0.81        -0.90     0.33 -0.09
## sodio           -0.10     0.43    0.93   0.81  1.00        -0.62     0.67 -0.07
## carboidratos    -0.59    -0.85   -0.64  -0.90 -0.62         1.00    -0.02  0.08
## calorias        -0.76     0.07    0.76   0.33  0.67        -0.02     1.00  0.02
## fds             -0.06    -0.05   -0.04  -0.09 -0.07         0.08     0.02  1.00
## borda            0.08     0.07    0.01   0.05  0.03        -0.08    -0.04  0.02
##              borda
## agua_100      0.08
## proteina      0.07
## gordura       0.01
## cinzas        0.05
## sodio         0.03
## carboidratos -0.08
## calorias     -0.04
## fds           0.02
## borda         1.00
## 
## n= 300 
## 
## 
## P
##              agua_100 proteina gordura cinzas sodio  carboidratos calorias
## agua_100              0.0000   0.0029  0.0000 0.0769 0.0000       0.0000  
## proteina     0.0000            0.0000  0.0000 0.0000 0.0000       0.2250  
## gordura      0.0029   0.0000           0.0000 0.0000 0.0000       0.0000  
## cinzas       0.0000   0.0000   0.0000         0.0000 0.0000       0.0000  
## sodio        0.0769   0.0000   0.0000  0.0000        0.0000       0.0000  
## carboidratos 0.0000   0.0000   0.0000  0.0000 0.0000              0.6854  
## calorias     0.0000   0.2250   0.0000  0.0000 0.0000 0.6854               
## fds          0.2728   0.4210   0.4408  0.1042 0.2442 0.1697       0.7862  
## borda        0.1563   0.2446   0.8210  0.3546 0.6387 0.1799       0.4428  
##              fds    borda 
## agua_100     0.2728 0.1563
## proteina     0.4210 0.2446
## gordura      0.4408 0.8210
## cinzas       0.1042 0.3546
## sodio        0.2442 0.6387
## carboidratos 0.1697 0.1799
## calorias     0.7862 0.4428
## fds                 0.7343
## borda        0.7343

Identificamos que algumas variáveis possuem correlação:

library(corrplot)
res1=cor(tabela_pizza)
corrplot(res1)

library(corrplot)
res1=cor(tabela_pizza[,-c(8,9)])
corrplot(res1,method = "square")

corrplot(res1,method = "square",order = "hclust",tl.srt = 45,pch=9,tl.col = "black",type="upper")

chart.Correlation(tabela_pizza[,-c(8,9)],histogram = T,pch=19)

col<- colorRampPalette(c("blue", "white", "red"))(20)

heatmap(x=res1,col=col,symm = T)

corrplot(res1,method = "pie")

corrplot(res1,method = "number")

if(!require(GGally)){install.packages("GGally")}
## Loading required package: GGally
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
ggpairs(tabela_pizza, lower = list(continuous = "smooth"))

ggpairs(pizza, columns = 3:7, ggplot2::aes(colour=marca))

ggcorr(res1, label=T)

Vamos analizar gordura e sodio.

cor.test(gordura,sodio)
## 
##  Pearson's product-moment correlation
## 
## data:  gordura and sodio
## t = 44.875, df = 298, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9170029 0.9465271
## sample estimates:
##       cor 
## 0.9333252

Como o valor de \(r\) é significativo, vamos construir o modelo linear.

plot(gordura,sodio)

modelo_pizza=lm(sodio~gordura,data=tabela_pizza)
anova(modelo_pizza)
summary(modelo_pizza)
## 
## Call:
## lm(formula = sodio ~ gordura, data = tabela_pizza)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.36215 -0.08424  0.01552  0.08489  0.44868 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.1096655  0.0189876  -5.776 1.93e-08 ***
## gordura      0.0385113  0.0008582  44.875  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1332 on 298 degrees of freedom
## Multiple R-squared:  0.8711, Adjusted R-squared:  0.8707 
## F-statistic:  2014 on 1 and 298 DF,  p-value: < 2.2e-16

Assim nosso modelo é expresso pela fórmula \[sodio_{prev}=-0.1096+0.0385*gordura + \epsilon\], com um \[r^2 = 0.8711\]

Então podemos ver que o modelo explica 87,11% da relação entre a gordura e o sódio na fabricação da pizza.

par(mfrow=c(2,2))
plot(modelo_pizza)

Vamos apresentar o pacotes easyanls que permita a construção de modelo lineares e não lineares entre duas variáveis numéricas.

if(!require(easynls)){install.packages("easynls")}
## Loading required package: easynls
## construindo a base de dados
dados_nls=data.frame(gordura,sodio)
nlsfit(dados_nls, model = 1,start = c(1,1))
## $Model
## [1] "y~a+b*x"
## 
## $Parameters
##                              sodio
## coefficient a          -0.10966552
## coefficient b           0.03851129
## p-value t.test for a    0.00000002
## p-value t.test for b    0.00000000
## r-squared               0.87110000
## adjusted r-squared      0.87070000
## AIC                  -354.21620831
## BIC                  -343.10486088

podendo também fazer o gráfico do modelo

nlsplot(dados_nls, model = 1, position = 6,
        xlab = "gordura",ylab="sódio")

Você pode defir o tipo de modelo a ser construido.

1 = "y~a+b*x" linear

2 = “y~a+bx+cx^2” quadratic

3 = “y ~ a + b * (x - c) * (x <= c)” linear plateau

4 = “y ~ (a + b * x + c * I(x^2)) * (x <= -0.5 * b/c) + (a + I(-b^2/(4 * c))) * (x > -0.5 * b/c)” quadratic plateau

5 = "ifelse(x>=d,(a-cd)+(b+c)x, a+b*x)" two linear

6 = “y~aexp(bx)” exponential

7 = "y~a(1+b(exp(-c*x)))^-1" logistic

8 = "y~a(1-b(exp(-c*x)))^3" van bertalanffy

9 = "y~a(1-b(exp(-c*x)))" brody

10 = "y~aexp(-bexp(-c*x)" gompertz

11 = "y~(ax^b)exp(-c*x)" lactation curve

12 = “y ~ a + b * (1 - exp(-c * x))” ruminal degradation curve

13 = “y~(a/(1+exp(2-4c(x-e))))+(b/(1+exp(2-4d(x-e))))” logistic bi-compartmental

nlsplot(dados_nls, model = 2, position = 6,
        xlab = "gordura",ylab="sódio")