Correlação e Regressão

Correlação

Reginaldo K Fukuchi

No dia a dia, termos como dependência, associação e correlação são comumente usadas como sinônimos. Entretanto, tecnicamnente, associação e correlação deveriam ser usados para expressar situações distintas. Associação seria uma relação mais geral entre variáveis, ou seja, uma variável fornece informação para outra variável. Por outro lado correlação é mais específica e significa que as variáveis estão correlacionadas apresentando padrões crescentes e decrescentes [1].

Correlation is a type of association and measures increasing or decreasing trends quantified using correlation coefficients. (a) Scatter plots of associated (but not correlated), non-associated and correlated variables. In the lower association example, variance in y is increasing with x. (b) The Pearson correlation coefficient (r, black) measures linear trends, and the Spearman correlation coefficient (s, red) measures increasing or decreasing trends. (c) Very different data sets may have similar r values. Descriptors such as curvature or the presence of outliers can be more specific [4].

A análise de correlação é um método estatístico para examinar a relação entre variáveis. Quando estudamos a relação entre apenas duas variáveis dizemos que se trata de uma correlação simples. Por exemplo, um pesquisador pode estar interessado em investigar se níveis de força muscular diminuem com a idade. Para isto ele coleta dados de força muscular em indivíduos de diversas faixas etárias na tentativa de responder esta questão. Nesta caso, a variável independente (ou explanatória) é a idade, enquanto a variável dependente é o nível de força muscular. Em alguns casos existem mais de uma variável independente e, portanto, usamos a chamada correlação múltipla. Quando estudamos o comportamente de duas variáveis a correlação pode ser positiva, negativa ou sem correlação. O exemplo abaixo mostra dados de revendedoras de automóveis onde o interesse é estudar se existe relação entre o número de veículos disponíveis e o lucro da empresa [2].

# Exemple of positive correlation.
carNumber <- c(63,29,20.8,19.1,13.4,8.5) # independent variable
revenue <- c(7,3.9,2.1,2.8,1.4,1.5) # dependent variable

# Showing results
plot(carNumber,revenue, type = "p", main = "Number of Cars vs. Revenue", xlab="Cars(in 10,000)", ylab="Revenue (billions)")

Podemos perceber que a relação é positiva, ou seja, existe a tendência de elevar o lucro a medida que aumenta o número de carros. Por outro lado podemos ter a situação inversa, ou seja, a cada incremento na variável independente existe a tendência de redução na variável dependente. O exemplo abaixo mostra a relação entre a ausência nas aulas e o desempenho na prova de uma amostra de alunos.

# Example of negative correlation.
numMisses <- c(6,2,15,9,12,5,8)
grade <- c(82,86,43,74,58,90,78)

# Showing results
plot(numMisses,grade, type = "p", main = "Frequency vs. Grade", xlab="Number of Absences", ylab="Final grade")

Embora qualitativamente os dois casos ilustram bem os tipos de correlação (positiva e negativa), como podemos saber quão relevante são estas correlações. Para isto precisamos obter um valor que dá a medida de importância da relação linear entre duas varáveis conhecida como coeficiente de correlação.

Coeficiente de correlação (Pearson’s product-moment correlation)

O coeficiente de correlação é uma medida da relação entre as variáveis e o seu valor pode se encontrar entre -1 e 1. Em outras palavras, o coeficiente de correlação, pode ser definido como a quantificação do grau (e direção) pelo qual duas varáveis aleatórias estão relacionadas, dado que esta relação seja linear. Quanto mais próximo dos extremos maior é a força da relação seja ela positiva (1) ou negativa (-1). Em geral, este coeficiente é simbolizado pela letra \(r\) para amostras e a letra grega \(\rho\) para a população. O coeficiente de correlação é calculado pela equação.

\[ r = \frac{n(\sum xy)-(\sum x)(\sum y)}{\sqrt{[n(\sum x^2)-(\sum x)^2][n(\sum y^2)-(\sum y^2)]}} \]

Onde n é o número observações, x é a variável independente e y a variável dependente.

Se rearranjarmos a equação acima obtemos:

\[ r = \frac{1}{(n-1)} \sum_{i=1}^{n} (\frac{x_i-\bar{x}}{S_x}) (\frac{y_i-\bar{y}}{S_y}) \]

Note que o r agora consiste simplesmente na média dos produtos dos scores \(Z\) para x e y.

Usando os dados do exemplo anterior das revendedoras de veículos podemos calcular o \(r\) no R.

x <- c(63,29,20.8,19.1,13.4,8.5) # independent variable
y <- c(7,3.9,2.1,2.8,1.4,1.5) # dependent variable
n <- 6 # sample size

# Correlation coefficient
r <- (n*(sum(x*y)) - sum(x)*sum(y))/sqrt((n*sum(x^2)-sum(x)^2)*(n*sum(y^2)-sum(y)^2))
print(paste0("The coefficient of correlation is ",toString(round(r,3))))

## [1] "The coefficient of correlation is 0.982"

Podemos usar a função cor.test da biblioteca “psych” no R para fazer isto e obter o mesmo resultado.

library(psych)
# Pairwise Pearson's product-moment correlation
r <- cor.test(x,y)
print(paste0("The coefficient of correlation is ",toString(round(r$estimate,3))))

## [1] "The coefficient of correlation is 0.982"

Portanto, podemos observar que existe uma relação substancial entre o número de veículos e o lucro das revendedoras. Entretanto, não podemos ainda dizer se esta relação é significativa pois a análise de correlação foi realizada em uma amostra da população (que é pequena por sinal) e os resultados, portanto, podem ter sido obtidos por acaso. Então, para determinar se o resultado é significativo ou não, usamos o teste de hipóteses. No caso da análise de correlação, estabelecemos duas hipóteses da seguinte forma:

• \(H_0\) de que não existe correlação entre variáveis independente e dependente
• \(H_1\) de que existe correlação significativa entre as variáveis

Portanto, se a \(H_0\) for rejeitada com um nível de significância, podemos dizer que o r é diferente de zero e que existe uma correlação significativa entre as variáveis. Entretanto para que o coeficiente de correlação obtido de uma amostra possa ser considerado representativo de uma população, é necessário que as seguintes condições sejam satisfeitas como descrito em [1]

1. As variáveis são linearmente relacionadas
2. As variáveis são aleatórias
3. As variáveis seguem uma distribuição normal bivariada.

Como as variáveis foram obtidas do mesmo local, elas não podem ser consideradas independentes. Portanto, para cada valor correspondente de x a distribuição de y será normal e vice-versa.

Teste de significância de r

Para testar a significância do r de uma amostra (que é tipicamente o caso) podemos usar o teste t como mostra a fórmula abaixo.

\[ t = r \sqrt{\frac{n-2}{1-r^2}} \]

Onde n - 2 são os graus de liberdade.

Usando o exemplo anterior, vamos examinar se a correlação entre número de veículos e o lucro das revendedoras é significativa. Como hábito, deveríamos sempre começar declarando as hipóteses a priori.

• \(H_0\): O tamanho da frota de carros não está relacionada com o lucro da revendedora.
• \(H_1\): O tamanho da frota de carros está relacionada com o lucro da revendedora.

Considerando um alpha de 0,05 e o r de 0,982, temos.

# Data
alpha = 0.05
r = 0.982
n = 6

# Applying in the equation
t = r*sqrt((n-2)/(1-r^2))
print(paste0("The t test is ",toString(round(t,3))))

## [1] "The t test is 10.398"

Precisamos, agora, encontrar o valor crítico para decidir se devemos rejeitar a $ H_0 $.

# Finding critical value given alpha level 0.05
cv <- qt(alpha/2,n-2, lower.tail=FALSE)
print(paste0("The critical value for the test is ",toString(round(cv,3))))

## [1] "The critical value for the test is 2.776"

Como o valor crítico é menor do que o valor do teste (10.398 vs. 2.776), temos subsídios para rejeitar \(H_0\) e concluir que existe sim uma correlação significativa entre o número de veículos e o lucro das revendedoras.

Lembre que também podemos tomar a decisão baseado no valor de p como segue.

ptest <- pt(t,n-2,lower.tail = FALSE)
print(paste0("The p value for the test is ",toString(ptest*2)))

## [1] "The p value for the test is 0.000483084000000001"

Portanto, chegamos ao mesmo resultado pois a probabilidade de rejeitar \(H_0\) quando esta é verdadeira é muito menor do que o alpha. Concluímos que o tamanho da frota de veículos das revendedores tem correlação positiva significativa com o lucro das mesmas. Na verdade, a função cor.test faz tudo isto como podemos observar no seu output.

# Pearson's product moment correlation coefficient
r <- cor.test(x,y, method = "pearson", alternative = "two.sided")
r

## 
##  Pearson's product-moment correlation
## 
## data:  x and y
## t = 10.392, df = 4, p-value = 0.0004842
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8392388 0.9981103
## sample estimates:
##       cor 
## 0.9819798

Regressão

“The statistician knows…that in nature there never was a normal distribution, there never was a straight line, yet with normal and linear assumptions, known to be false, he can often derive results which match, to a useful approximation, those found in the real world.” Cited in [1].

Ao observar que existe uma relação linear entre duas variáveis (por meio da análise de correlação), um pesquisador pode se interessar agora em examinar a tendência dos dados e a capacidade de fazer predições a partir destes dados. Em um sentido amplo, a análise de regressão tenta ajustar uma reta que represente o comportamento dos dados da melhor maneira possível.

Fonte: (https://madhureshkumar.wordpress.com/2015/07/21/regression-analysis-linear-regression-sse-assumption-of-linear-regression-error-term-and-best-fit-line/)

Se as variáveis não estão relacionadas, não existe razão para conduzir uma análise de regressão. A figura abaixo ilustra os diferentes casos quando duas variáveis são estudadas.

Figura.(a) If the properties of Y do not change with X, there is no association. (b) Association is possible without regression. Here E(Y|X) is constant, but the variance of Y increases with X. (c) Linear regression E(Y|X) = β0+ β1X. (d) Nonlinear regression E(Y|X) = exp(β0+ β1X). Source [1]

A reta que melhor se ajusta aos dados é aquela que minimiza a soma das distâncias verticais ao quadrado de cada ponto a esta reta. Portanto, esta reta é determinada por um método conhecido como mínimos quadrados (least squares method)

A reta que melhor se ajusta aos dados pode ser representada matematicamente pela equação abaixo:

\[ \hat{Y} = \hat{\beta{_0}} + \hat{\beta{_1}} X \]

Onde, \(\hat{\beta{_0}}\) e \(\hat{\beta{_1}}\) é o intercepto (coeficiente linear) e inclinação (coeficiente angular) da reta, respectivamente. Portanto, o objetivo da regressão é encontrar \(\hat{\beta{_0}}\) e \(\hat{\beta{_1}}\) de modo que a soma quadrática dos desvios de cada observação a reta ajustada seja a mínima possível.

Figura. Desvio entre o valor observado e o valor predito pelo modelo de regressão. Fonte [2].

A reta que melhor se ajusta aos dados pode portanto ser calculada da seguinte forma [3].

Os desvios entre o valor observado e o predito pelo modelo representa os resíduos. Portanto, o modelo de regressão, incluindo os resíduos seria.

\[ \hat{Y} = \hat{\beta{_0}} + \hat{\beta{_1}} X + \epsilon \]

Onde, \(\epsilon = observado - predito\).

Então, a qualidade do ajuste do modelo tem relação com estes desvios entre o valor predito e observado.

Variação do modelo de regressão

Quanto menor a diferença entre os valores observados e os preditos melhor é o modelo e, tipicamente, a correlação entre as variáveis também será boa ou próxima de 1 ou -1 [2]. Portanto, alguns tipos de variações no modelo de regressão precisam ser compreendidas. A variação total do modelo de regressão consiste na diferença entre a variação explicada e não explicada pelo modelo de acordo com a equação abaixo e ilustrado na Figura.

\[\sum{(y-\bar{y})^2} = \sum{(y'-\bar{y})^2} + \sum{(y-y')^2}\]

Ou seja \[SS_T = SS_R + SS_E\]

Onde SS é a soma dos quadrados. Figura. Ilustração das variâncias do modelo de regressão. Fonte [2].

A reta horizontal na média de y representa a condição onde não existe correlação entre x e y. Portanto, o \(SS_T\) indica o desvio total de y pois o x não contribuiu para este modelo. Por outro lado, quando x participa do modelo para predizer y, o desvio são os resíduos. Portanto, para examinar a qualidade do modelo no ajuste aos dados utilizamos a razão entre a varição explicada pelo modelo pela variação total que também é conhecido como coeficiente de determinação (\(R^2\)) do modelo de regressão linear.

\[R^2 = \frac{SS_T-SS_E}{SS_T} = \frac{SS_R}{SS_T} = \frac{Var. Explicada}{Var. Total}\]

Em outras palavras, o \(R^2\) representa a variação na variável dependente explicada pelo modelo (linha de regressão) e pela variável independente.

Vamos calcular estas variações usando os seguintes dados hipotéticos do exemplo 10-3 [2].

x <- c(1,2,3,4,5) # independent variable
y <- c(10,8,12,16,20) # dependent variable
n <- 5 # sample size
ybar <- mean(y) # average of y

# Visualizing data
plot(y~x, main = "Scatterplot", ylab = "Dependent variable (y)", xlab = "Independent variable (x)")

# Calculating regression coefficients
coeff <- lm(y~x) # linear fit to get coefficients of regression model

# Getting coefficients and calculating the predicted values
yhat <- coeff$coefficients[1] + coeff$coefficients[2]*x

abline(coeff$coefficients[1],coeff$coefficients[2], col = "red") # plot the regression line

A equação da reta que se ajusta aos pontos é \(y'=\) 4.8 \(+\) 2.8 \(x\)

Calculando a variação total (\(SS_T\)), variação explicada (\(SS_R\)) e não explicada (\(SS_E\)) pelo modelo.

# Total variation
SSt <- sum((y-ybar)^2)
print(paste0("The total sum of squares is ",toString(SSt)))

## [1] "The total sum of squares is 92.8"

# Variation explained
SSr <- sum((yhat-ybar)^2)
print(paste0("The explained variation is ",toString(round(SSr,2))))

## [1] "The explained variation is 78.4"

# Variation unexplained
SSe <- sum((y-yhat)^2)
print(paste0("The unexplained variation is ",toString(round(SSe,2))))

## [1] "The unexplained variation is 14.4"

Agora podemos calcular então o coeficiente de determinação \(R^2\)

rSquared <- (SSt - SSe)/SSt
print(paste0("The coeff. of determination is ",toString(round(rSquared,2))))

## [1] "The coeff. of determination is 0.84"

Homocedasticidade

Se lembrarmos que as diferenças entre os valores preditos e observados contituem os resíduos do modelo, podemos examinar a natureza da distribuição destes resíduos em torno da média que é zero. Se a distribuição dos resíduos é normal podemos sugerir que uma relação entre as variáveis (independente e dependente) existe e que o desvio padrão da variável dependente em cada valor correspondente da variável independente é o mesmo (Figura). Se os dados possuem estas características podemos dizer que eles atendem a condição de homocedasticidade ou homogeneidade das variâncias. Figura. In a linear regression relationship, the response variable has a distribution for each value of the independent variable. (a) At each height, weight is distributed normally with s.d. s = 3. (b) Linear regression of n = 3 weight measurements for each height. The mean weight varies as m(Height) = 2 × Height/3 – 45 (black line) and is estimated by a regression line (blue line) with 95% confidence interval (blue band). The 95% prediction interval (gray band) is the region in which 95% of the population is predicted to lie for each fixed height. Fonte [1].

Para verificar a distribuição dos resíduos podemos visualizá-los por meio de gráficos de distribuição. Vamos calcular os resíduos no R.

# Calculating residuals
res <- y-yhat

# Visualing residuals.
plot(res~x, main = "Residuals plot", ylab = "Residuals", xlab = "Independent variable")
abline(0,0, col = "red") # horizontal line representing mean residuals (equal zero)

Podemos observar que os dados se dispersam de forma normal em torno da média e, portanto, indica que a homocedasticidade foi satisfeita. Existem testes estatísticos para verificar esta propriedade como o Levene’s test e o Bartlet’s test mas os mesmos não serão apresentados aqui. A figura abaixo mostra algumas condições onde a distribuição dos resíduos não é tão favorável para fazer predições.

Erro padrão e intervalo de predição

Nos exemplos anteriores, tratamos da predição de um valor da variável dependente de acordo com modelo de regressão mas é uma predição pontual. Entretanto, assim como vimos anteriormente, também é possível prever um intervalo no qual o valor da variável dependente de uma porcentagem da população dado que a variável dependente está fixa. Este intervalo é conhecido como intervalo de predição. Por exemplo, na relação entre a altura e o peso corporal, a faixa cinza representa o intervalo de predição e a faixa azul o intervalo de confiança. Perceba que o intervalo de confiança é mais estreito pois consiste na variação da estimativa da média do peso corporal com 95% de confiança. Para o cálculo do intervalo de predição, é necessário obter o erro padrão da estimativa da predição como mostrado abaixo.

\[S_{est} = \sqrt{\frac{\sum{(y-y')^2}}{n-2}}\]

Note que o cálculo acima é similar ao do desvio padrão mas não calculamos o desvio da média mas sim o desvio do modelo de regressão. Portanto, quanto mais próximo os valores preditos estão dos observados menor o erro padrão da predição. O intervalo de predição pode então ser calculado da seguinte forma.

\[ y'-t_{ \alpha /2} S_{est} \sqrt{1+\frac{1}{n}+\frac{n(x-\bar{X})^2}{n \sum{x^2}-(\sum{x})^2}} <y< y'+t_{ \alpha /2} S_{est} \sqrt{1+\frac{1}{n}+\frac{n(x-\bar{X})^2}{n \sum{x^2}-(\sum{x})^2}}\]

Vamos calcular o erro padrão da média e o intervalo de predição usando o exemplo do livro do Bluman.

Example 10-12. A researcher collects the following data and determines that there is a significant relationship between the age of a copy machine and its monthly maintenance cost.

# Data set
age <- c(1,2,3,4,4,6)
cost<- c(62,78,70,90,93,103)
n <- 6

# Calculating regression coefficients
coef <- lm(cost~age)

Portanto, a equação da reta que se ajusta aos dados é: \(y'=\) 55.57 \(+\) 8.13 \(x\)

# Calculating the predicted values
yhat <- coef$coefficients[1] + coef$coefficients[2]*age

# Standard error of the estimate calculation
Sest <- sqrt(sum((cost-yhat)^2)/(n-2))
print(paste0("The std. error of the estimate is ",toString(round(Sest,2))))

## [1] "The std. error of the estimate is 6.51"

Vamos calcular agora o intervalo de predição.

#Calculating the critical value based on t distribution
df <- n-2
alpha <- 0.05
cv <- qt(alpha/2,df,lower.tail=FALSE)
# Lower and upper limit of prediction interval given the age
lowerLim <- yhat-cv*Sest*sqrt(1+(1/n)+(n*(age-mean(age))^2)/(n*sum(age^2)-sum(age)^2))
upperLim <- yhat+cv*Sest*sqrt(1+(1/n)+(n*(age-mean(age))^2)/(n*sum(age^2)-sum(age)^2))

E o intervalo de confiança.

#95% Confidence interval calculation

# Lower and upper limit of prediction interval given the age
lowerBound <- yhat-cv*Sest*sqrt((1/n)+(n*(age-mean(age))^2)/(n*sum(age^2)-sum(age)^2))
upperBound <- yhat+cv*Sest*sqrt((1/n)+(n*(age-mean(age))^2)/(n*sum(age^2)-sum(age)^2))

Vamos visualizar os resultados

# Visualizing data distribution
plot(cost~age, main = "Age vs Maintainance cost", ylab = "Maintainance cost ($)", xlab = "Age (years)")

# Visualizing the regression line
abline(coef$coefficients[1],coef$coefficients[2], col = "red") 

# Visualizing lower and upper prediction interval only for the experimental observations
lines(age,lowerLim,type = "l", col= "blue")
lines(age,upperLim,type = "l", col= "blue")

# Confidence interval
lines(age,lowerBound,type = "l", col= "black")
lines(age,upperBound,type = "l", col= "black")

Perceba que o intervalo de predição é mais largo do que o intervalo de confiança. Isso se deve pelo fato que o intervalo de predição objetiva prever observações futuras então além do erro associado ao modelo de regressão ele também precisa incluir o erro associado a predições de observações futuras.

References

1. Altman & Krzywinski (2015).Points of Significance: Simple linear regression. Nature Methods 12, 999–1000.
2. Bluman, Allan G. Elementary statistics : a step by step approach / Allan Bluman. — 8th ed.
3. Montgomery (2013) Applied Statistics and Probability for Engineers. John Wiley & Sons.
4. Altman & Krzywinski (2015). Points of Significance: Association, correlation and causation. Nature Methods 12, 899–900