Regressão Linear

Agora que já aprendemos ‘tudo’ sobre ANOVA, vamos aprender um pouco sobre um outro tipo de modelo linear que é usado para analisar a relação entre variáveis contínuas: a Regressão Linear. A regressão linear é usada quando: (1) se quer saber se uma variável contínua está associada a outra variável contínua; (2) quando se quer medir a força da associação (r2); ou (3) se quer a equação que descreve a relação para poder usá-la na predição de valores que não são conhecidos.

A regressão linear (simples ou múltipla) é feita com a função lm(), a exemplo do ANOVA - que também é um modelo linear. Relembrando, essa função requer uma fórmula. No caso da regressão linear simples a fórmula assume a forma DV ~ IV, que podemos ler como ‘DV como função de IV’ ou ‘DV predita por IV’, ‘DV modelada por IV’, etc. Lembre-se que DV (ou variável resposta) deve vir antes de ‘~’ e IV ou variáveis explanatórias depois. É super simples, vamos tentar usando a base de dados ‘cats’ do pacote ‘MASS’, que contém informação sobre algumas características de gatos domésticos.

require(MASS)
## Loading required package: MASS
data(cats)
str(cats)
## 'data.frame':    144 obs. of  3 variables:
##  $ Sex: Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Bwt: num  2 2 2 2.1 2.1 2.1 2.1 2.1 2.1 2.1 ...
##  $ Hwt: num  7 7.4 9.5 7.2 7.3 7.6 8.1 8.2 8.3 8.5 ...
summary(cats)
##  Sex         Bwt             Hwt       
##  F:47   Min.   :2.000   Min.   : 6.30  
##  M:97   1st Qu.:2.300   1st Qu.: 8.95  
##         Median :2.700   Median :10.10  
##         Mean   :2.724   Mean   :10.63  
##         3rd Qu.:3.025   3rd Qu.:12.12  
##         Max.   :3.900   Max.   :20.50

‘Bwt’ é a massa corpórea em kilogramas, ‘Hwt’ é a massa do coração em gramas, em machos e fêmeas (‘Sex’). Vamos checar o comportamento dos dados usando um scatterplot.

attach(cats)
plot(Bwt, Hwt)
title(main="Massa do Coração (g) vs. Massa corpórea (kg)\nde Gatos Domésticos")

A função plot( ) retorna um scatterplot sempre que alimentada com duas variáveis numéricas. A primeira variável listada é plotada no eixo horizontal. Também é possível usar a fórmula DV ~ IV, e o resultado é o mesmo.

plot(Hwt ~ Bwt, main="Massa do Coração (g) vs. Massa corpórea (kg)\nde Gatos Domésticos")

O gráfico mostra uma relação que parece ser forte e razoavelmente linear entre a massa do coração e a massa corpórea. Vamos testar usando lm()

lm(Hwt ~ Bwt)
## 
## Call:
## lm(formula = Hwt ~ Bwt)
## 
## Coefficients:
## (Intercept)          Bwt  
##     -0.3567       4.0341

Daí podemos tirar a equação da regressão: Hwt = 4.0341 (Bwt) - 0.3567. Mas para ter acesso a mais informações, é importante armazenar o output na forma de um objeto e extrair dele outras informações com outras funções.

fit <- lm(Hwt ~ Bwt)
fit
## 
## Call:
## lm(formula = Hwt ~ Bwt)
## 
## Coefficients:
## (Intercept)          Bwt  
##     -0.3567       4.0341
# reparem que essa função extrai muito mais informação
# incluindo o teste estatístico
summary(fit) 
## 
## Call:
## lm(formula = Hwt ~ Bwt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.5694 -0.9634 -0.0921  1.0426  5.1238 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -0.3567     0.6923  -0.515    0.607    
## Bwt           4.0341     0.2503  16.119   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.452 on 142 degrees of freedom
## Multiple R-squared:  0.6466, Adjusted R-squared:  0.6441 
## F-statistic: 259.8 on 1 and 142 DF,  p-value: < 2.2e-16
anova(fit)
## Analysis of Variance Table
## 
## Response: Hwt
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## Bwt         1 548.09  548.09  259.83 < 2.2e-16 ***
## Residuals 142 299.53    2.11                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nós podemos plotar a reta de regressão com o output do objeto para nos ajudar na interpretação. Para publicar nós vamos usar o ggplot2, que permite muito mais customização e fica bem mais bonito!

plot(Hwt ~ Bwt, main="Gatinhos")
abline(fit, col="red")

Há vários diagnósticos de regressão disponíveis, o mais simples dele pode ser feito usando a função plot()

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

O primeiro comando que usamos ajusta os parâmetros gráficos do dispositivo, particionando a janela gráfica em 4 partes (2 linhas e 2 colunas). Dessa forma, quatro gráficos podem ser plotados em apenas uma janela. A janela gráfica vai reter essa estrutura até que você a feche (no RSTUDIO seria o equivalente a usar o botão da vassoura na janela de plots), retornando então ao default. Outra opção é usar par(mfrow=c(1,1)).

O primeiro plot de diagnósticos é um standard residual plot mostrando os resíduos contra os valores ajustados. Pontos que apresentam uma tendência outlier são identificados. O que esperamos é que não apareça nenhuma tendência nos pontos deste plot, caso contrário o modelo linear pode não ser o mais adequado para estudar essa relação. O segundo plot é um normal quantile plot, que mostra se a distribuição dos resíduos segue ou não uma normal (esperamos que sim!). O último plot mostra residuals vs. leverage, cujos pontos identificados podem estar influenciando demais a tendência da regressão. Um deles é o caso 144.

cats[144,]
##     Sex Bwt  Hwt
## 144   M 3.9 20.5
fit$fitted[144]
##      144 
## 15.37618
fit$residuals[144]
##      144 
## 5.123818

Ele era um gato gordo (3.9kg) e tinha um coração enorme (20.5g)! Os maiores de todo o dataset. O valor observado de ‘Hwt’ tem uma diferença (resíduo) de 5.1g em relação ao ajustado (15.4g). O erro padrão residual (vejo o output do modelo acima) foi 1.452. Convertendo para resíduo normalizado temos 5.124/1.452=3.53, e isso é muito! Uma medida comum de influência seria Cook’s Distance, vamos tentar.

par(mfrow=c(1,1))                 
plot(cooks.distance(fit))

Qual o valor mais discrepante? E agora o que fazer? Uma coisa que podemos fazer é ajustar um novo modelo sem o gato 144 e comparar os coeficientes.

fit.sem.144 <- lm(Hwt ~ Bwt, subset=(Hwt<20.5))
fit.sem.144
## 
## Call:
## lm(formula = Hwt ~ Bwt, subset = (Hwt < 20.5))
## 
## Coefficients:
## (Intercept)          Bwt  
##       0.118        3.846

Há também outras opções, como ajustar um modelo que é mais robusto tendo em vista que o dataset apresenta outliers. Existe uma função no pacote MASS chamada rlm() que é interessante para estes casos.

ANCOVA

A análise de covariância (ANCOVA) é usada para comparar duas ou mais retas de regressão, testando o efeito de um fator (variável categórica) em uma variável dependente (y) enquanto controla o efeito de uma co-variável contínua (x).

As retas são comparadas estudando a interação de uma variável categórica (por exemplo tratamentos, sexo) com a variável independente. Vamos tentar com os dados que já estávamos trabalhando. A pergunta que queremos testar é se machos e fêmeas são diferentes quanto ao tamanho do coração. Não podemos usar uma ANOVA pois o tamanho do coração depende do tamanho do corpo, como vimos no exemplo anterior.

fit2 <- aov(Hwt ~ Bwt * Sex)
summary(fit2)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## Bwt           1  548.1   548.1 263.645 <2e-16 ***
## Sex           1    0.2     0.2   0.074 0.7853    
## Bwt:Sex       1    8.3     8.3   4.008 0.0472 *  
## Residuals   140  291.0     2.1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Na ANCOVA, a primeira hipótese nula que testamos é a de que os coeficientes das retas de regressão (neste caso a reta das fêmeas e a dos machos) são iguais, ou seja as retas são paralelas. Isso é feito analisando a interação entre o fator e a covariável. Se a interação é significativamente diferente de zero, os efeitos da covariável na resposta dependem do nível do fator. Ou seja, as retas de regressão têm diferentes coeficientes de inclinação. Neste caso, vemos que o efeito do aumento no tamanho do corpo na massa do coração é distinto entre machos e fêmeas. Com a interação significativa, não podemos seguir com o teste da segunda hipótese nula da ANCOVA, a de que os interceptos das retas são iguais (os grupos não são diferentes para tratamento, sexo…). Temos que fazer análises separadas para cada sexo e podemos ver o resultado obtido de forma gráfica.

machos <- subset(cats, Sex=="M")
femeas <- cats[cats$Sex=='F',]

fitm <- lm(Hwt~Bwt, data=machos)
summary(fitm)
## 
## Call:
## lm(formula = Hwt ~ Bwt, data = machos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7728 -1.0478 -0.2976  0.9835  4.8646 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.1841     0.9983  -1.186    0.239    
## Bwt           4.3127     0.3399  12.688   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.557 on 95 degrees of freedom
## Multiple R-squared:  0.6289, Adjusted R-squared:  0.625 
## F-statistic:   161 on 1 and 95 DF,  p-value: < 2.2e-16
fitf <- lm(Hwt~Bwt, data=femeas)
summary(fitf)
## 
## Call:
## lm(formula = Hwt ~ Bwt, data = femeas)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.00871 -0.68599 -0.04506  0.79583  2.21858 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.9813     1.4855   2.007 0.050785 .  
## Bwt           2.6364     0.6254   4.215 0.000119 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.162 on 45 degrees of freedom
## Multiple R-squared:  0.2831, Adjusted R-squared:  0.2671 
## F-statistic: 17.77 on 1 and 45 DF,  p-value: 0.0001186
plot(Hwt~Bwt, type='n')
points(machos$Bwt,machos$Hwt, pch=20, col="blue")
points(femeas$Bwt,femeas$Hwt, pch=1, col='red')
abline(fitm, lty=2, col='blue')
abline(fitf, lty=1, col='red')
legend("bottomright", c("Machos","Fêmeas"), lty=c(2,1), 
       pch=c(20,1), col=c('blue', 'red') )

Com um outro dataset, vamos testar se o número de sementes que uma determinada planta produz é alterado pela qualidade do solo, dada pelo tipo de fertilizante usado no experimento. Como sabemos que a produtividade dessa planta é uma função da chuva, usamos a quantidade de água oferida a planta durante o experimento como covariável.

sementes <- read.csv('http://renatabrandt.github.io/EBC2015/data/sementes.csv')
sem.fit1 = lm(Sementes ~ Chuva*Fertilizante, data = sementes)
anova(sem.fit1)
## Analysis of Variance Table
## 
## Response: Sementes
##                      Df Sum Sq Mean Sq    F value Pr(>F)    
## Chuva                 1  36068   36068 34423.6954 <2e-16 ***
## Fertilizante          2  24952   12476 11907.1841 <2e-16 ***
## Chuva:Fertilizante    2      1       1     0.7147 0.4895    
## Residuals          1494   1565       1                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

O resultado nos mostra que a quantidade de sementes produzidas depende da quantidade de água oferecida para a planta. Além disso, que essa relação não difere entre os tratamentos de fertilizantes. Podemos prosseguir com a análise para testar nossa hipótese principal retirando o termo de interação.

sem.fit2 = lm(Sementes ~ Chuva + Fertilizante, data = sementes)
anova(sem.fit2)
## Analysis of Variance Table
## 
## Response: Sementes
##                Df Sum Sq Mean Sq F value    Pr(>F)    
## Chuva           1  36068   36068   34437 < 2.2e-16 ***
## Fertilizante    2  24952   12476   11912 < 2.2e-16 ***
## Residuals    1496   1567       1                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(sem.fit1,sem.fit2)
## Analysis of Variance Table
## 
## Model 1: Sementes ~ Chuva * Fertilizante
## Model 2: Sementes ~ Chuva + Fertilizante
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1   1494 1565.4                           
## 2   1496 1566.9 -2   -1.4978 0.7147 0.4895

O teste nos mostra que remover a interação não afeta o ajuste do modelo significantemente. Logo, ficamos com o modelo mais parsimonioso. Biologicamente então temos que a quantidade de sementes produzidas pelas nossas plantas aumenta com a quantidade de água fornecida para as plantas (proxy para chuva), e que os tratamentos com diferentes fertilizantes aumentam a produção de sementes. Mas queremos ver visualmente certo?

a <- subset(sementes, Fertilizante=="F1")
b <- subset(sementes, Fertilizante=="F2")
c <- subset(sementes, Fertilizante=="F3")

fita <- lm(Sementes ~ Chuva, data=a)
summary(fita)
## 
## Call:
## lm(formula = Sementes ~ Chuva, data = a)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7399 -0.6119  0.0064  0.7138  3.2488 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 104.11322    1.36555   76.24   <2e-16 ***
## Chuva         4.96353    0.04545  109.21   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9954 on 498 degrees of freedom
## Multiple R-squared:  0.9599, Adjusted R-squared:  0.9598 
## F-statistic: 1.193e+04 on 1 and 498 DF,  p-value: < 2.2e-16
fitb <- lm(Sementes ~ Chuva, data=b)
summary(fitb)
## 
## Call:
## lm(formula = Sementes ~ Chuva, data = b)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.08953 -0.65306  0.01996  0.67215  2.91962 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 97.56039    1.38130   70.63   <2e-16 ***
## Chuva        5.01556    0.04597  109.10   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.007 on 498 degrees of freedom
## Multiple R-squared:  0.9598, Adjusted R-squared:  0.9598 
## F-statistic: 1.19e+04 on 1 and 498 DF,  p-value: < 2.2e-16
fitc <- lm(Sementes ~ Chuva, data=c)
summary(fitc)
## 
## Call:
## lm(formula = Sementes ~ Chuva, data = c)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.89939 -0.69463  0.01832  0.72268  2.84367 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 106.78492    1.46411   72.94   <2e-16 ***
## Chuva         5.04106    0.04873  103.45   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.067 on 498 degrees of freedom
## Multiple R-squared:  0.9555, Adjusted R-squared:  0.9554 
## F-statistic: 1.07e+04 on 1 and 498 DF,  p-value: < 2.2e-16
attach(sementes)
plot(Sementes ~ Chuva, type='n', ylab = 'Produção de Sementes',
     xlab = 'Água fornecida durante o experimento', ylim = c(225,280))
points(a$Chuva,a$Sementes, pch=21, bg="blue")
points(b$Chuva,b$Sementes, pch=21, bg='red')
points(c$Chuva,c$Sementes, pch=21, bg="purple")

abline(fita, lty=1, col='blue')
abline(fitb, lty=1, col='red')
abline(fitc, lty=1, col='purple')

legend("bottomright", c("F1","F2", "F3"), lty= 1, 
       pch=21, col=c('blue', 'red', 'purple') )

Uma terceira possibilidade para um outcome do ANCOVA seria o efeito do tratamento e interação com a covariável não significativos. Se o efeito da covariável é significativo, significa que há uma única reta de regressão. No caso das sementes, a produção aumentaria com a quantidade de água fornecida durante o experimento mas não haveria diferença entre os tratamentos de fertilizantes.

’Otimo!!! Ja possuímos os resultados dos testes e os graficos que podemos salvar em .pdf. Mas como construir uma tabela adequada para apresentar em um artigo? Com o R em mãos você não vai ficar copiando e colando cada resultado no excel, vai?

Nós temos a solução para o seu problema! Stargazer package!!!! / / O pacote stargazer foi idealizado exatamente com o intuito de construir tabelas para publicação sem muito sofrimento.

Ele pode retornar as tabelas em formato .txt, .csv, entre outros. Entretanto nós recomendamos salvar em formato .html. Como exemplo, vamos usar os resultados de fita, fitb & fitc relativos a produtivadade das plantas.

install.packages(“stargazer”)

require(stargazer)
## Loading required package: stargazer
## 
## Please cite as: 
## 
##  Hlavac, Marek (2015). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2. http://CRAN.R-project.org/package=stargazer
stargazer (fita, fitb, fitc, title="O STARGAZER É DEMAIS", 
              align=T, type="html", style="all", digits.extra=4, single.row=T)
## 
## <table style="text-align:center"><caption><strong>O STARGAZER É DEMAIS</strong></caption>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="3"><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="3" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td colspan="3">Sementes</td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Chuva</td><td>4.964<sup>***</sup> (0.045)</td><td>5.016<sup>***</sup> (0.046)</td><td>5.041<sup>***</sup> (0.049)</td></tr>
## <tr><td style="text-align:left"></td><td>t = 109.210</td><td>t = 109.097</td><td>t = 103.450</td></tr>
## <tr><td style="text-align:left"></td><td>p = 0.000</td><td>p = 0.000</td><td>p = 0.000</td></tr>
## <tr><td style="text-align:left">Constant</td><td>104.113<sup>***</sup> (1.366)</td><td>97.560<sup>***</sup> (1.381)</td><td>106.785<sup>***</sup> (1.464)</td></tr>
## <tr><td style="text-align:left"></td><td>t = 76.242</td><td>t = 70.629</td><td>t = 72.935</td></tr>
## <tr><td style="text-align:left"></td><td>p = 0.000</td><td>p = 0.000</td><td>p = 0.000</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>500</td><td>500</td><td>500</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.960</td><td>0.960</td><td>0.956</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.960</td><td>0.960</td><td>0.955</td></tr>
## <tr><td style="text-align:left">Residual Std. Error (df = 498)</td><td>0.995</td><td>1.007</td><td>1.067</td></tr>
## <tr><td style="text-align:left">F Statistic (df = 1; 498)</td><td>11,926.850<sup>***</sup> (p = 0.000)</td><td>11,902.170<sup>***</sup> (p = 0.000)</td><td>10,701.830<sup>***</sup> (p = 0.000)</td></tr>
## <tr><td colspan="4" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="3" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>

Apos “rodar” a funcao devemos copiar o texto gerado no console abaixo e colar no bloco de notas, salvando o arquivo como .html. Pronto! Se clicarmos no arquivo ele abrirá em formato de tabela. Você pode ver a tabela que construímos aqui. Outra opção é “Abrir com -> Microsoft Word” e sua tabela estará lá! O stargazer pode ser aplicado pra mutos tipos de dados e modelos, desde tabelas descritivas até testes estatísticos filogenéticos.

That’s all folks!