Correlação e Regressão

A relação entre variáveis quantitativas contínuas, como as medidas tomadas entre os caranguejos do dataset crabs, pode ser avaliado através do cálculo da correlação linear, cujos índices mais conhecidos são o de Pearson e o de Spearman, os quais são calculados com o mesmo comando cor. Neste caso, vamos ver se existe uma relação entre o tamanho do lobo frontal (FL) e a largura da carapaça (CW) dos caranguejos.

library(MASS)
attach(crabs)
cor(crabs$FL, crabs$CW)
## [1] 0.9649558

O primeiro índice de correlação é o de Pearson que é o default do comando cor. Se quiser usar outro índice de correlação, este é informado através do comando: method=“spearman”.

cor(crabs$FL, crabs$CW, method="spearman")
## [1] 0.9654331

Observe os resultados de ambos os índices de correlação. Note quão alta é associação entre estas duas variáveis através do plot.

plot(FL~CW, data=crabs)

Esta relação entre as duas variáveis pode ser descrita através de uma regressão linear. A regressão é um caso de Modelo Linear (“Linear Model”), cujo comando é lm. Assim:

minha.regressao<-lm(FL~CW, data=crabs)
summary(minha.regressao)
## 
## Call:
## lm(formula = FL ~ CW, data = crabs)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.17541 -0.73802 -0.02847  0.77317  1.76698 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.019233   0.308461  -0.062     0.95    
## CW           0.428462   0.008281  51.743   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9195 on 198 degrees of freedom
## Multiple R-squared:  0.9311, Adjusted R-squared:  0.9308 
## F-statistic:  2677 on 1 and 198 DF,  p-value: < 2.2e-16
anova(minha.regressao)
## Analysis of Variance Table
## 
## Response: FL
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## CW          1 2263.83 2263.83  2677.4 < 2.2e-16 ***
## Residuals 198  167.42    0.85                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
minha.regressao$coefficients
## (Intercept)          CW 
## -0.01923255  0.42846208

No comando summary você tem os dados da regressão como o intercepto (-0.02) e a declividade (0.43), o valor de R-quadrado=0.93, além da significância da regressão: p<2.2*10^-16

Uma tabela mais sintética é obtida com o comando anova indicando a significância da regressão. Utilizando-se do comando abline podemos agora acrescentar, no plot, o modelo de regressão, ou seja, a linha azul que melhor explica a relação entre FL e CW.

plot(FL~CW, data=crabs)
abline(minha.regressao, lwd=2, col='blue')

O comando abline é aplicado ao modelo gerado que é a equação da reta, os parâmetros lwd e col indicam a espessura e cor da linha, estão aqui apenas para mostrar que você pode alterar os parâmetros gráficos. Caso você opte por não incluir estes comandos, ficando apenas com abline(minha.regressão), vai ser plotada uma linha mais fina e preta, que é o padrão (default) do plot. Para saber qual o padrão do plot você pode usar o help:

?plot

Como nós vimos, apesar de termos uma reta de ajuste (regressão) parece haver dois padrões principais pela distribuição dos pontos. Como são dados de comprimento frontal e largura de carapaça de duas espeçies diferentes (sp O e sp B), talvez os dois padrões sejam destas espécies. Para isso o ideal seria fazer uma regressão de cada um deles separadamente e plotar as duas no mesmo gráfico. Como?

Primeiro separando os dados de cada espécie , isto é fazendo um subset dos dados:

spB<-subset(crabs, sp=="B")
spO<-subset(crabs, sp=="O")

Depois, fazendo duas regressões, reg1 e reg2, uma para cada espécie:

reg1<-lm(FL~CW, data=spB)
reg2<-lm(FL~CW, data=spO)

E, por útlimo, plotando novamente a figura e colocando as duas retas de ajuste das regressões separadamente:

plot(FL~CW, data=crabs)
abline(reg1, lwd=2, col='blue')
abline(reg2, lwd=2, col='red')

Notamos, portanto que realmente aquele padrão que, aparentemente, era duplo se refere mesmo as duas espécies cujas regressões são diferentes.

Análise de Variância e Teste t

Como visto na descrição dos dados (comando ?crabs), FL e CW são medidas obtidas de duas espécies de caranguejos identificados em machos e fêmeas. Vimos também que há um padrão diferente entre estas duas espécies mas, e cada variável individualmente? Por exemplo, será que a largura da carapaça (variável CW) é diferente nas duas espécies (variável sp)?

plot(CW~sp, data=crabs)

A figura acima é um boxplot, veja na literatura como interpretar. Note que o mesmo comando plot fornece figuras diferentes. Isto aconteceu, pois agora estamos lidando com uma variável numérica (CW) e uma variável categórica, ou fator (sp). No plot anterior da regressão lidávamos com duas variáveis quantitativas (numéricas), por isso fizemos uma regressão. Agora que temos fatores categóricos vamos fazer uma Análise de Variância usando os mesmos comandos, o próprio R vai reconhecer que sp é um fator, já que ele é composto por letras que representam as duas espécies (e não números). É recomendável quando você construir sua planilha que os códigos dos fatores sejam na forma de letras ou nomes e não de números, caso opte por números você terá que “avisar” o R que aquela variável é um fator e não uma variável quantitativa.

minha.anova<-lm(CW~sp, data=crabs)
anova(minha.anova)
## Analysis of Variance Table
## 
## Response: CW
##            Df  Sum Sq Mean Sq F value   Pr(>F)   
## sp          1   576.3  576.30  9.7069 0.002108 **
## Residuals 198 11755.3   59.37                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

O resultado mostra que existe sim uma diferença significativa entre as espécies quanto à largura da carapaça. Se você pedir o summary ou minha.anova$coefficientes , verá o mesmo resultado e mais algumas informações:

summary(minha.anova)
## 
## Call:
## lm(formula = CW ~ sp, data = crabs)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.5120  -5.7383   0.7855   5.5830  19.8830 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  34.7170     0.7705  45.057  < 2e-16 ***
## spO           3.3950     1.0897   3.116  0.00211 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.705 on 198 degrees of freedom
## Multiple R-squared:  0.04673,    Adjusted R-squared:  0.04192 
## F-statistic: 9.707 on 1 and 198 DF,  p-value: 0.002108
minha.anova$coefficients
## (Intercept)         spO 
##      34.717       3.395

Neste resultado o intercepto (34.717) indica o valor médio do comprimento da carapaça para a espécie B como vemos:

mediaSpB<-mean(crabs$CW[1:100])
mediaSpB
## [1] 34.717
mediaSpO<-mean(crabs$CW[101:200])
mediaSpO
## [1] 38.112

Note que, para efetuar a média da espécie “B” (mediaSpB) eu selecionei apenas os 100 primeiros valores usando a notação [1:100], enquanto que, para efetuar a média da espécie “O” (mediaSpB) eu usei as 100 restantes [101:200].

O valor médio da espécie “B” é exatamente o mesmo (34.717) do intercepto da saída dos coeficientes da “minha.anova”, sendo que a média da espécie “O” (os valores de 101 a 200) foi maior. Se subtrairmos uma média da outra, teremos:

mediaSpO-mediaSpB
## [1] 3.395

O valor da diferença obtido é o próprio coeficiente (slope) da análise (3.395). Este indica, portanto, o efeito do fator sp (espécie) na estimativa da média das carapaças. Bom, agora você sabe interpretar os dados. Note que você só pode fazer esta interpretação pois trabalhou com os dados brutos. Muitas vezes, determinadas análises requerem algum tipo de transformação dos dados para que as premissas da análise estatística sejam válidas. Outra opção, que veremos mais tarde, é o uso de outras análises mais adequadas para aquele tipo de distribuição dos dados.
Caso seja necessária uma transformação, você pode fazer direto na própria linha de comando, sem ter que alterar totalmente a sua planilha original. Assim, a mesma Análise de Variância efetuada aqui pode ser feita com os dados de largura da carapaça (CW) em escala logarítmica ou extraindo sua raiz quadrada. Isto é feito colocando na frente da variável que está entre parênteses a função do R para estas operações: log(CW), para logaritmo ou sqrt(CW), para extrair a raiz quadrada.

minha.anova.log<-lm(log(CW)~sp, data=crabs)
minha.anova.raiz<-lm(sqrt(CW)~sp, data=crabs)

Neste caso fizemos uma análise de variância onde testamos se a largura da carapaça (CW) era diferente entre as duas espécies (O e B). Como a variável sp (espécies) só tem dois níveis, isto é as espécies B e O, isso pode ser feito também através de um teste t:

t.test(CW~sp, data=crabs)
## 
##  Welch Two Sample t-test
## 
## data:  CW by sp
## t = -3.1156, df = 197.65, p-value = 0.002109
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.54389 -1.24611
## sample estimates:
## mean in group B mean in group O 
##          34.717          38.112

Reparem que a saída do teste t é mais sintética e os valores são praticamente idênticos ao da análise de variância, a diferença se deve apenas a algum arredondamento. Ou seja, para uma variável quantitativa com dois níveis a anova e o teste t são equivalentes.

Análise de variância fatorial (ou de duas vias)

A análise de variância fatorial permite trabalhar com dois (idealmente) ou mais fatores (raramente) na mesma análise. Neste caso, vimos como a largura da carapaça dos caranguejos (CW) variou entre as espécies, mas temos também diferenças entre sexos (variável: sex). Vamos ver como são estas diferenças entre espécies (sp) e sexo (sex) ao mesmo tempo?

minha.anova2<-lm(CW~sp*sex, data=crabs)
anova(minha.anova2)
## Analysis of Variance Table
## 
## Response: CW
##            Df  Sum Sq Mean Sq F value   Pr(>F)   
## sp          1   576.3  576.30 10.0567 0.001762 **
## sex         1    68.3   68.33  1.1924 0.276196   
## sp:sex      1   455.1  455.11  7.9419 0.005325 **
## Residuals 196 11231.8   57.31                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Note o **”*“** entre os dois fatores (sp e sex) indicando a interação. Se você preferir fazer uma anova de duas vias sem interação deve usar o símbolo **”+“**.

O resultado indica que, além de uma variação no comprimento entre espécies diferentes, há também uma interação relacionada ao sexo. Para entender esta interação você pode fazer o seguinte plot:

 interaction.plot(crabs$sp,crabs$sex,crabs$CW)

Observe que além das espécies serem diferentes quanto a largura da carapaça, na espécie “B” os machos são maiores, enquanto que na espécie “O” eles são menores que as fêmeas, por isso há uma interação significativa. Esta interação significativa aparece na tabela acima da anova na linha sp:sex.

Faça você a mesma análise, mas sem interação (substituindo o **”*“** por **”+“**) e compare as tabelas da anova no resultado

minha.anova2<-lm(CW~sp+sex, data=crabs)
anova(minha.anova2)

As análises efetuadas através de modelos lineares, no caso tanto a Regressão quanto a ANOVA, devem ser efetuadas quando os dados estão adequados às premissas destas análises. Você deve ver em algum livro básico de estatística ou bioestatística quais são estas premissas (por ex: distribuição normal, homogeneidade das variâncias). Entretanto é possível, após realizar a análise, fazer um diagnóstico da qualidade através da plotagem de um modelo. Para maiores detalhes clique no link Diagnostico

Análise de Covariância

Quando fizemos a regressão lá em cima, nós plotamos duas regressões, uma para cada grupo, mas quando fizemos o modelo linear (lm) testamos apenas se havia uma relação entre o tamanho do lobo frontal (FL) e a largura da carapaça (CW) dos caranguejos e notamos, graficamente que esta era diferente ente as duas espécies. O teste para mostrar se esta diferença, que parece clara, é significativa é o que denominamos de Análise de Covariância onde comparamos se o padrõa de dependência entre CW e FL é diferente entre as espécies.

No caso temos uma variável dita dependente e contínua, FL e uma indepedente ou explicativa também contínua, que é CW e uma variável independente explicativa que é categórica sp. O modelo linear da ANCOVA é, portanto, uma mistura da análise de regressão (variável explicativa contínua) e a análise de variância (variável explicativa categórica) que vimos acima.

minha.covariancia<-lm(FL~CW*sp, data=crabs)
anova(minha.covariancia)
## Analysis of Variance Table
## 
## Response: FL
##            Df  Sum Sq Mean Sq   F value    Pr(>F)    
## CW          1 2263.83 2263.83 16856.242 < 2.2e-16 ***
## sp          1  134.17  134.17   999.015 < 2.2e-16 ***
## CW:sp       1    6.92    6.92    51.549 1.414e-11 ***
## Residuals 196   26.32    0.13                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Neste caso vamos ver linha a linha a saída: (1) na primeira linha vemos que há uma dependência de FL por CW, como já visto na análise de regressão, (2) na segunda linha vemos o efeito das espécies em FL ou se o Lobo Frontal (FL) é diferente nas duas espécies para uma mesma largura (CW). Neste caso existe uma diferença no valor do intercepto das duas espécies, ou seja o valor delas no eixo Y quando o valor de X é igual a 0. (3) Na terceira linha temos a interação entre CW e sp, esta interação indica que a inclinação (slope) de cada uma das espécies é diferente (o valor de p tambem é significativo).

plot(FL~CW, data=crabs)
abline(reg1, lwd=2, col='blue')
abline(reg2, lwd=2, col='red')

Portanto, tanto o intercepto quanto a inclinação das regressões são signficativos. Este é o objetivo principal da análise de covariância ver se o padrão de dependência observado entre as duas variáveis contínuas é o mesmo para duas categorias, no caso as duas espécies.

Acima vimos casos onde se avalia se avalia a dependência entre uma variável contínua por outras variáve contínua (regressão), uma variável contínua por uma variável categórica (ANOVA, Teste t) ou duas ou mais variáveis categóricas (ANOVA Fatorial). Mas, e quando nenhuma das variáveis é contínua, todas são categóricas? Vamos ver abaixo.

Qui-quadrado

O teste de qui-quadrado procura avaliar se duas variáveis categóricas são independentes, ou seja, se uma delas afeta a outra. O exemplo abaixo consta do pacote MASS e envolve os hábitos de fumar de estudantes e o nível de atividade física dos mesmos.

library (MASS)
head(survey$Smoke)
## [1] Never Regul Occas Never Never Never
## Levels: Heavy Never Occas Regul
head(survey$Exer)
## [1] Some None None None Some Some
## Levels: Freq None Some

Acima estão mostrados apenas os cinco primeiros casos (comando head), se quiser ver os demais basta digitar o nome da variável. Note que as variáveis acima são categóricas - Smoke indicando o hábito de fumar em níveis decrescentes (níveis: Heavy, Occas, Regul, Never) e hábito de exercício, Exer (Freq, Some, None) Agora vamos transformar estes dados categóricos em uma tabela de contingência usando o comando table. A tabela de contingênica é uma tabela de combinação de duas variáveis categóricas a qual fornece quantos casos correspondem a cada combinação. No nosso caso vai fornecer a informação de quantos casos há em cada combinação de Smoke e Exer.

tabela.qui<-table(survey$Smoke, survey$Exer)
tabela.qui
##        
##         Freq None Some
##   Heavy    7    1    3
##   Never   87   18   84
##   Occas   12    3    4
##   Regul    9    1    7

Usando esta tabela vamos fazer o teste do qui-quadrado:

chisq.test(tabela.qui)
## Warning in chisq.test(tabela.qui): Chi-squared approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  tabela.qui
## X-squared = 5.4885, df = 6, p-value = 0.4828

O resultado não foi significativo (ver o valor de p), indicando que não houve uma relação, ao menos entre os estudantes avaliados, entre prática esportiva e hábito de fumar. Note que existe um aviso (Warning) dizendo que os resultados podem não estar corretos. Isto se deve ao fato de que os valores de algumas células da tabela são muito baixos. Por exemplo, as combinações Heavy e None e Regul e None. Nestes casos existem outras análises mais adequadas, como o Teste Exato de Fisher:

fisher.test(tabela.qui)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  tabela.qui
## p-value = 0.4138
## alternative hypothesis: two.sided

Veja que desta vez não houve aviso, indicando que o teste é adequado e, novamente, que não existe relação entre prática esportiva e hábito de fumar (ao menos neste grupo estudado).

Continuação: Gráficos