---
title: "CE083 - Estatística Computacional I"
author: "Seu nome"
date: "17 de dezembro, 2015"
output: html_document
---

# Sobre R Markdown

- Veja o material disponível em [Rmarkdown-rautu](http://git.leg.ufpr.br/fernandomayer/Rmarkdown-rautu)
- Veja apresentação do [R Users Group](http://www.leg.ufpr.br/~fernandomayer/useR/Rmarkdown-rautu.html)
- Veja o vídeo da apresentação no [Youtube](https://www.youtube.com/watch?v=qFh4XtckP2k)

# Regressão linear

Importa a base de dados `crabs.csv`

```{r}
url <- "http://www.leg.ufpr.br/~fernandomayer/dados/crabs.csv"
```

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

```{r}

```

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
  
- Exemplo:
	- Sabendo a altura ($X$) podemos determinar o peso ($Y$) de uma
      pessoa, se conhecemos os parâmetros do modelo anterior
	- Sabendo a medida CL ($X$) de um caranguejo, podemos determinar
      (**predizer**) qual seria sua medida CW ($Y$)

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:
  
- $Y$ é a variável **resposta** (ou dependente)
- $X$ é a variável **explicativa** (ou independente)
- $\beta_0$ é o **intercepto** da reta (valor de $Y$ quando
  $X = 0$)
- $\beta_1$ é o **coeficiente angular** da reta
  (**efeito** de $X$ sobre $Y$)
- $e_i \sim \text{N}(0, \sigma^2)$ é o **erro**, ou
    **desvio**, ou **resíduo**

**Interpretação dos parâmetros:**

- $\beta_0$ representa o ponto onde a reta corta o eixo $Y$ (na maioria
das vezes não possui interpretação prática)
- $\beta_1$ representa a variabilidade em $Y$ causada pelo aumento de
uma unidade em $X$. Além disso,
	- $\beta_1 > 0$ mostra que com o aumento de $X$, também há um
      aumento em $Y$
	- $\beta_1 = 0$ mostra que **não há efeito** de $X$ sobre $Y$
	- $\beta_1 < 0$ mostra que com a aumento de $X$, há uma diminuição
      em $Y$
   
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

```{r}

```

Agora calculamos $\hat{\beta_1}$ com

```{r}

```

E $\hat{\beta_0}$ é calculado com

```{r}

```

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

```{r}

```

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

```{r}

```

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

```{r}

```

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()`
  
```{r}

```

#### 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

```{r, eval=FALSE}
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:

```{r}

```

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

```{r}

```

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

```{r}

```

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

```{r}

```

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()`

```{r}

```

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

```{r}

```

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

```{r}

```

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()`

```{r}

```

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

```{r}

```

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

```{r}

```

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,

```{r, eval=FALSE}
sqrt(mod.an$Sum[2]/mod.an$Df[2])
```

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

```{r, eval=FALSE}
mod.an$F[1]
```

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

```{r, eval=FALSE}
mod.sum$coef[2,3]^2
```

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.

```{r}

```

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

```{r}

```

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

```{r}

```

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

```{r}

```

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()`

```{r}

```

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

```{r, fig.show='hold'}
## Deixa a janela gráfica com 2 linhas e 2 colunas

```

## 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)
