Sobre R Markdown

Regressão linear

Importa a base de dados crabs.csv

url <- "http://www.leg.ufpr.br/~fernandomayer/dados/crabs.csv"
dados <- read.table(url, header = TRUE, sep = ";", dec = ",")
str(dados)
## 'data.frame':    156 obs. of  7 variables:
##  $ especie: Factor w/ 2 levels "azul","laranja": 1 1 1 1 1 1 1 1 1 1 ...
##  $ sexo   : Factor w/ 2 levels "F","M": 2 2 2 2 2 2 2 2 2 2 ...
##  $ FL     : num  8.1 8.8 9.2 9.6 10.8 11.6 11.8 12.3 12.6 12.8 ...
##  $ RW     : num  6.7 7.7 7.8 7.9 9 9.1 10.5 11 10 10.9 ...
##  $ CL     : num  16.1 18.1 19 20.1 23 24.5 25.2 26.8 27.7 27.4 ...
##  $ CW     : num  19 20.8 22.4 23.1 26.5 28.4 29.3 31.5 31.7 31.5 ...
##  $ BD     : num  7 7.4 7.7 8.2 9.8 10.4 10.3 11.4 11.4 11 ...

Vamos analisar a relação que existe entre CL e CW

plot(CW ~ CL, data = dados)

Um modelo linear entre duas variáveis \(X\) e \(Y\), é definido matematicamente como uma equação com dois parâmetros desconhecidos,

\[ Y = \beta_0 + \beta_1 X \]

A análise de regressão é a técnica estatística que analisa as relações existentes entre uma única variável dependente, e uma ou mais variáveis independentes O objetivo é estudar as relações entre as variáveis, a partir de um modelo matemático, permitindo estimar o valor de uma variável a partir da outra

O problema da análise de regressão consiste em definir a forma de relação existente entre as variáveis. Por exemplo, podemos ter as seguintes relações

\[ \begin{align} Y &= \beta_0 + \beta_1 X &\qquad \text{linear} \\ Y &= \beta_0 X^{\beta_1} &\qquad \text{potência} \\ Y &= \beta_0 e^{\beta_1 X} &\qquad \text{exponencial} \\ Y &= \beta_0 + \beta_1 X + \beta_2 X^2 &\qquad \text{polinomial} \\ \end{align} \]

Em todos os casos, a variável dependente é \(Y\), aquela que será predita a partir da relação e da variável independente \(X\).

Em uma análise de regressão linear consideraremos apenas as variáveis que possuem uma relação linear entre si. Uma análise de regressão linear múltipla pode associar \(k\) variáveis independentes (\(X\)) para “explicar” uma única variável dependente (\(Y\)),

\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_k X_k + e \]

Uma análise de regressão linear simples associa uma única variável independente (\(X\)) com uma variável dependente (\(Y\)),

\[ Y = \beta_0 + \beta_1 X + e \]

Assim, dados \(n\) pares de valores, \((X_1, Y_1), (X_2, Y_2), \ldots, (X_n, Y_n)\), se for admitido que \(Y\) é função linear de \(X\), pode-se estabelecer uma regressão linear simples, cujo modelo estatístico é

\[ Y_i = \beta_0 + \beta_1 X_i + e_i, \quad i = 1, 2, \ldots, n \]

onde:

Interpretação dos parâmetros:

O problema agora consiste em estimar os parâmetros \(\beta_0\) e \(\beta_1\).

Estimação dos parâmetros

Como através de uma amostra obtemos uma estimativa da verdadeira equação de regressão, denominamos

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

ou seja, \(\hat{Y}_i\) é o valor estimado de \(Y_i\), através das estimativas de \(\beta_0\) e \(\beta_1\), que chamaremos de \(\hat{\beta}_0\) e \(\hat{\beta}_1\). Para cada valor de \(Y_i\), temos um valor \(\hat{Y}_i\) estimado pela equação de regressão,

\[ Y_i = \hat{Y}_i + e_i \]

Portanto, o erro (ou desvio) de cada observação em relação ao modelo adotado será

\[ \begin{align} e_i &= Y_i - \hat{Y}_i \\ e_i &= Y_i - (\beta_0 + \beta_1 X_i) \end{align} \]

Devemos então adotar um modelo cujos parâmetros \(\beta_0\) e \(\beta_1\), tornem esse diferença a menor possível. Isso equivale a minimizar a soma de quadrados dos resíduos (\(SQR\)), ou do erro,

\[ SQR = \sum_{i=1}^{n} [Y_i - (\beta_0 + \beta_1 X_i)]^2 \]

O método de minimizar a soma de quadrados dos resíduos é denominado de método dos mínimos quadrados. Para se encontrar o ponto mínimo de uma função, temos que obter as derivadas parciais em relação a cada parâmetro,

\[ \begin{align} \frac{\partial SQR}{\partial \beta_0} &= 2 \sum_{i=1}^{n} [Y_i - \beta_0 - \beta_1 X_i] (-1) \\ \frac{\partial SQR}{\partial \beta_1} &= 2 \sum_{i=1}^{n} [Y_i - \beta_0 - \beta_1 X_i] (-X_i) \end{align} \]

e igualar os resultados a zero

\[ \hat{\beta}_0 = \frac{\partial SQR}{\partial \beta_0} = 0 \qquad \text{e} \qquad \hat{\beta}_1 = \frac{\partial SQR}{\partial \beta_1} = 0 \]

Dessa forma, chegamos às estimativas de mínimos quadrados para os parâmetros \(\beta_0\) e \(\beta_1\):

\[ \begin{align} \hat{\beta}_1 &= \frac{\sum_{i=1}^{n} X_iY_i - \frac{\sum_{i=1}^{n} X_i \sum_{i=1}^{n} Y_i}{n}}{\sum_{i=1}^{n}X_i^2 - \frac{(\sum_{i=1}^{n} X_i)^2}{n}} \\ & \\ \hat{\beta_0} &= \bar{Y} - \hat{\beta}_1 \bar{X} \end{align} \]

onde

\[ \begin{align} \bar{Y} = \frac{1}{n} \sum_{i=1}^{n} Y_i \qquad \text{e} \qquad \bar{X} = \frac{1}{n} \sum_{i=1}^{n} X_i \end{align} \]

Portanto, podemos usar estas soluções para calcular os parâmetros \(\hat{\beta_0}\) e \(\hat{\beta_1}\) na relação linear entre CL e CW, ou seja,

\[ CW = \hat{\beta}_0 + \hat{\beta}_1 CL + e \]

Estimando os parâmetros à mão

Como vimos pelas soluções acima, primeiro calculamos \(\hat{\beta_1}\), e depois \(\hat{\beta_0}\). Para facilitar as contas, vamos criar objetos X e Y com as colunas CL e CW respectivamente, e n que é o tamanho da amostra

X <- dados$CL
Y <- dados$CW
n <- length(X)

Agora calculamos \(\hat{\beta_1}\) com

beta1 <- (sum(X * Y) - (sum(X) * sum(Y))/n)/(sum(X^2) - (sum(X)^2/n))
beta1
## [1] 1.097451

E \(\hat{\beta_0}\) é calculado com

beta0 <- mean(Y) - beta1 * mean(X)
beta0
## [1] 1.18695

Uma parte importante em uma análise de regressão linear é a verificação dos resíduos do modelo, ou seja, os desvios de cada valor observado \(Y\) em relação aos valores preditos pelo modelo, \(\hat{Y}\).

Os valores preditos (ou estimados) de \(Y\) pelo modelo podem ser calculados usando os coeficientes estimados

Y.pred <- beta0 + beta1 * X

Portanto, os resíduos podem ser calculados como a diferença entre cada valor observado, e cada valor predito de \(Y\), ou seja,

res <- Y - Y.pred

Como vimos que a suposição do modelo é de que os resíduos possuam uma distribuição normal com média 0 e variância constante, \(e_i \sim \text{N}(0, \sigma^2)\), então podemos verificar essa suposição fazendo um histograma destes resíduos

hist(res)

Note que por esta inspeção visual, a média dos resíduos parece estar próxima de 0, e todos os valores estão entre -2 e 2, o que indica que não há valores extremos em relaçao à normal padrão.

Observação: Existem outros métodos para analisar os resíduos e verificar as suposições do modelo, mas essa não é a intenção aqui.

Usando funções prontas para estimar os parâmetros

Para ajustar um modelo linear no R, usamos a função lm()

mod <- lm(CW ~ CL, data = dados)
mod
## 
## Call:
## lm(formula = CW ~ CL, data = dados)
## 
## Coefficients:
## (Intercept)           CL  
##       1.187        1.097

Sobre o uso de fórmulas no R

Entre os diversos usos de fórmulas, o mais importante deles é sem dúvida o fato que fórmulas são utilizadas na declaração de modelos estatísticos. Um aspecto particularmente importante da linguagem S, o portanto no programa R, é que adota-se uma abordagem unificada para modelagem, o que inclui a sintaxe para especificação de modelos. Variáveis respostas e covariáveis (variáveis explanatórias) são sempre especificadas de mesma forma básica, ou seja, na forma

resposta ~ covariavel

onde:

  • à esquerda indica-se a(s) variável(eis) resposta
  • o símbolo ~ significa “é modelada por”
  • à direita indica-se a(s) covariável(eis)

No restante deste texto vamos, por simplicidade, considerar que há apenas uma variável resposta que poderá ser explicada por uma ou mais covariáveis. Considere, o conjunto de dados do exemplo acima. Na sintaxe

lm(CW ~ CL, data = dados)

deve-se ler: CW é modelada por CL, através de um modelo linear lm(), o que implica no modelo acima.

Continuando…

Vimos acima que o objeto mod contém o modelo linear, e ao chamar esse objeto, vemos os coeficientes estimados, os mesmo que calculamos manualmente. Note a classe desse objeto:

class(mod)
## [1] "lm"

Existem uma série de funções genéricas capazes de extrair informações que estão nesse objeto. Objetos da classe lm possuem muito mais informações do que os coeficientes. Você pode conferir quais são estas informações com

names(mod)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

e para acessar algum destes componentes, fazemos como se fosse uma lista. Por exemplo, para acessar os coeficientes usamos

mod$coefficients
## (Intercept)          CL 
##    1.186950    1.097451

No entanto, também podemos usar a função genérica coef()

coef(mod)
## (Intercept)          CL 
##    1.186950    1.097451

Outras funções úteis nesse caso são:

  • fitted() para acessar os valores ajustados, ou preditos pelo modelo
  • residuals() para acessar os resíduos do modelo

Outra função genérica extremamente importante para modelos de qualquer classe é a função summary()

summary(mod)
## 
## Call:
## lm(formula = CW ~ CL, data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7762 -0.5699  0.1098  0.4629  1.8273 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.186950   0.285340    4.16 5.28e-05 ***
## CL          1.097451   0.008698  126.17  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7827 on 154 degrees of freedom
## Multiple R-squared:  0.9904, Adjusted R-squared:  0.9904 
## F-statistic: 1.592e+04 on 1 and 154 DF,  p-value: < 2.2e-16

Note que o resultado é um resumo completo sobre a distribuição dos resíduos do modelo, e uma tabela com os coeficientes estimados, incluindo um teste \(t\) para a hipótese nula de que os coeficientes são iguais a zero.

Veja também que podemos atribuir o sumário a um outro objeto, verificar sua classe e seus componentes

mod.sum <- summary(mod)
class(mod.sum)
## [1] "summary.lm"
names(mod.sum)
##  [1] "call"          "terms"         "residuals"     "coefficients" 
##  [5] "aliased"       "sigma"         "df"            "r.squared"    
##  [9] "adj.r.squared" "fstatistic"    "cov.unscaled"

Portanto, verificamos que outras informações estão disponíveis também para esse objeto, resultante do mesmo modelo. Por exemplo, para acessar o valor de \(R^2\) do modelo, usamos

mod.sum$r.squared
## [1] 0.9904188

Outra vantagem de se ajustar um modelo ocm a função lm() é que podemos também avaliar a significância geral do modelo através de uma Análise de Variância (ANOVA) para a regressão. Como vimos na estimação dos parâmetros, o objetivo é encontrar parâmetros que façam com que a soma de quadrados dos resíduos seja mínima. Podemos particionar a soma de quadrados da seguinte forma:

\[ SQTot = SQMod + SQRes \]

Portanto, se um modelo é bem ajustado, esperamos que a soma de quadrados do modelo \(SQMod\) seja grande, e a \(SQRes\) sejá mínima. Um quadro de ANOVA para o modelo irá testar, através de um teste F, se a soma de quadrados do modelo é significativamente diferente de zero. Para fazer essa ANOVA, usamos a função anova()

anova(mod)
## Analysis of Variance Table
## 
## Response: CW
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## CL          1 9752.6  9752.6   15919 < 2.2e-16 ***
## Residuals 154   94.3     0.6                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Com esse resultado, vemos que o modelo linear é significativo, ou seja, a variável explicativa CL é útil para explicar a relação com a variável resposta CW.

Podemos também associar a saída dessa ANOVA para um objeto, e verificar sua classe

mod.an <- anova(mod)
class(mod.an)
## [1] "anova"      "data.frame"

Veja que este objeto possui duas classes: anova e data.frame. Fundamentamente, a estrutura do objeto é a de um data.frame, portanto, podemos acessar essa tabela usando indexação conforme qualquer outro data frame. Veja os nomes das colunas e linhas desse objeto

names(mod.an)
## [1] "Df"      "Sum Sq"  "Mean Sq" "F value" "Pr(>F)"
row.names(mod.an)
## [1] "CL"        "Residuals"

Por exemplo, veja que o Residual standard error: 0.7827, como reportado no objeto mod.sum é o estimador do desvio-padrão residual \(\hat{\sigma}_{e} = \sqrt{\frac{\text{SQRes}}{n-2}}\), ou seja,

sqrt(mod.an$Sum[2]/mod.an$Df[2])
## [1] 0.7827079

e que F-statistic: 1.592e+04 (~15920) é o mesmo valor de

mod.an$F[1]
## [1] 15919.11

que testa a mesma hipótese da ANOVA. De fato, o valor de \(t^2\) para \(\beta_1\) no sumário do modelo é

mod.sum$coef[2,3]^2
## [1] 15919.11

Ainda aproveitando o objeto mod, podemos ajustar o modelo graficamente ao gráfico da relação entre CW e CL. Para isso, usamos uma função gráfica de baixo nível, abline(), utilizada para inserir linhas em gráficos. Esta função aplicada à um objeto da classe lm produz o seguinte resultado.

plot(CW ~ CL, data = dados)
abline(mod)

Apenas para ilustrar a interpretação dos parâmetros, vamos alterar a escala dos eixos X e Y para visualizar o intercepto do modelo

plot(CW ~ CL, data = dados, xlim = c(0,50), ylim = c(0,55))
abline(mod)

Assim como fizemos antes, vamos verificar a adequção do modelo através dos resíduos. Para isso, podemos extrair os resíduos diretamente do objeto mod com a função residuals(), e fazer o histograma diretamente

hist(residuals(mod))

Veja que os resíduos retornados pela função são exatamente os mesmos daqueles que calculamos anteriormente (objeto res)

all.equal(res, residuals(mod), check.names = FALSE)
## [1] TRUE

Isto é obviamente consequência do fato de que os valores preditos pelos coeficientes que calculamos à mão (objeto Y.pred) são os mesmos daqueles retornados pela função fitted()

all.equal(Y.pred, fitted(mod), check.names = FALSE)
## [1] TRUE

Uma forma mais completa de se analisar os resíduos de um modelo é através do método lm para a função genérica plot(). Ao utilizar essa função em um objeto da classe lm, por padrão são gerados quatro gráficos com diferentes informações sobre os resíduos

## Deixa a janela gráfica com 2 linhas e 2 colunas
par(mfrow = c(2, 2))
plot(mod)
par(mfrow = c(1, 1))

Exercícios

Com as colunas BD e CL do objeto dados:

  • Faça um gráfico da relação entre estas variáveis
  • Ajuste um modelo linear
    • Veja o sumário e o quadro de ANOVA
    • Ajuste a linha do modelo no gráfico
    • Verifique os resíduos
  • Qual sua conclusão?
    • Existe uma relação significativa? De que tipo (positiva, negativa)?
    • O modelo linear descreve bem a relação entre estas duas variáveis?
  • O modelos foi bem ajustado aos dados (observe os resíduos)