07-modelos_lineares.Rmd

# Modelos Lineares

A  classe de *Modelos Lineares* atende a uma ampla gama de problemas aplicados, apresentada em profundidade por [@neter2005applied]. Para uma introdução à classe de *Modelos Lineares Generalizados* recomenda-se [@mccullagh1989generalized].


## Correlação  {#correlacao}

<!-- ## Coeficiente de correlação (produto-momento) de Pearson -->
<!-- p. 1249 -->
<!-- ```{definition} -->
<!-- *Correlação* é o grau de associação linear entre duas variáves.  $\\$ -->
<!-- ``` -->

<!-- ## Coeficiente de correlação ordinal de Spearman -->
<!-- p. 1365 -->

<!-- ## Tau de Kendall??? -->
<!-- p. 1383 -->

<!-- ## Coeficiente de concordância de Kendall??? -->
<!-- p.1399 -->

## Regressão Linear Simples

### Modelo
O modelo de *Regressão Linear Simples* universal/populacional é construído, pela abordagem clássica, com todos os $N$ pares ordenados do universo, e pode ser descrito pela relação a seguir.

\begin{equation}
Y_i = \beta_0 + \beta_1 X_i + \varepsilon_i,
(\#eq:rls)
\end{equation}
onde $\varepsilon_i \sim \mathcal{N}(0,\sigma_{\varepsilon})$.

Na maioria dos casos práticos trabalha-se com amostras, sendo necessário estimar os valores de $\beta_0$ e $\beta_1$. O método dos *mínimos quadrados (ordinários)* é utilizado para calcular estas estimativas. O princípio do método é minimizar a soma de quadrado dos erros, i.e., 
\begin{equation}
minimizar \sum_{i=1}^n \varepsilon_{i}^{2}.
(\#eq:mqo)
\end{equation}

### Estimativas dos parâmetros
Basicamente utiliza-se $\varepsilon_{i} = Y_{i} - \beta_0 - \beta_1 X_i$ da Eq. \@ref(eq:rls) e deriva-se (parcialmente) em relação a $\beta_0$ e $\beta_1$, fazendo cada uma das derivadas parciais igual a zero. Para maiores detalhes recomenda-se [@degroot2012probability]. As estimativas por mínimos quadrados são enfim dadas por 
\begin{equation}
\hat{\beta}_1 = \frac{n \sum{x_i y_i} - \sum{x_i} \sum{y_i}}{n \sum{x_i^2} - (\sum{x_i})^2} 
(\#eq:beta1)
\end{equation}
e 
\begin{equation}
\hat{\beta}_0 = \bar{y} - \hat{\beta}_1 \bar{x}.
(\#eq:beta0)
\end{equation}

### Análise de diagnóstico

#### Teste para $\beta_0$ {-}
As hipóteses usuais do teste para $\beta_0$ são $H_0: \beta_0 = 0$ *vs* $H_1: \beta_0 \ne 0$. Sob $H_0$ 
\begin{equation}
T_0 = \frac{\hat{\beta}_0}{ep(\hat{\beta}_0)} \sim t_{n-2},
(\#eq:teste_b0)
\end{equation}
onde
\begin{equation}
ep(\hat{\beta}_0) = \sqrt{\hat{\sigma}^2 \left[ \frac{1}{n} + \frac{\bar{x}^2}{S_{xx}} \right] } = \sqrt{\frac{\sum_{i=i}^n (y_i - \hat{y}_i)^2}{n-2} \left[ \frac{1}{n} + \frac{\bar{x}^2}{\sum_{i=i}^n (x_i - \bar{x})^2} \right] }.
(\#eq:ep_b0)
\end{equation}

#### Teste para $\beta_1$ {-}

O teste para $\beta_1$ é fundamental na análise de diagnóstico. É com ele que decide-se a respeito da presença ou ausência de relação linear entre $X$ e $Y$. As hipóteses usuais do teste para $\beta_1$ são $H_0: \beta_1 = 0$ *vs* $H_1: \beta_1 \ne 0$. Sob $H_0$ 
\begin{equation}
T_1 = \frac{\hat{\beta}_1}{ep(\hat{\beta}_1)} \sim t_{n-2},
(\#eq:t1)
\end{equation}
 onde 
\begin{equation}
ep(\hat{\beta}_1) = \sqrt{\frac{\hat{\sigma}^2}{S_{xx}}} = \sqrt{\frac{\sum_{i=i}^n (y_i - \hat{y}_i)^2 / (n-2)}{\sum_{i=i}^n (x_i - \bar{x})^2} }.
(\#eq:ep_t1)
\end{equation}
 

### Modelo RPO
Um modelo na forma da Eq. \@ref(eq:rls0) é chamado *regressão pela origem* pelo fato de a reta ajustada passar pelo ponto $(0,0)$, a *origem* do plano cartesiano.
\begin{equation}
Y_i = \beta_1 X_i + \varepsilon_i,
(\#eq:rls0)
\end{equation}
onde $\varepsilon_i \sim \mathcal{N}(0,\sigma_{\varepsilon})$. A estimativa do parâmetro $\beta_1$ é dada por
\begin{equation}
\hat{\beta}_1 = \frac{\sum{x_i y_i} }{\sum{x_i^2} }.
(\#eq:beta1_rpo)
\end{equation}

#### Teste para $\beta_1$ do modelo RPO {-}

\begin{equation}
T_1^{RPO} = \frac{\hat{\beta}_1}{ep(\hat{\beta}_1)} \sim t_{n-1},
(\#eq:teste_beta1_rpo)
\end{equation}
onde
\begin{equation}
ep(\hat{\beta}_1) = \sqrt{\frac{\hat{\sigma}^2}{S_{xx}}} = \sqrt{\frac{\sum_{i=i}^n (y_i - \hat{y}_i)^2 / (n-1)}{\sum_{i=i}^n x_{i}^{2} } }.
(\#eq:ep_beta1_rpo)
\end{equation}


```{example}
O dono de um bar decidiu modelar o número de garrafas de bebidas vendidas em seu estabelecimento ($Y$) em função da temperatura máxima do dia ($X$).
```

```{r}
x <- read.table('http://www.filipezabala.com/data/drinks.txt', header = T, sep = '\t')
dim(x)  # dimensão de x
head(x) # cabeçalho de x
fit <- lm(gar ~ temp, data = x)       # regressão linear simples
summary(fit)
fit0 <- lm(gar ~ temp - 1, data = x)  # regressão pela origem
summary(fit0)
```


```{exercise}
Considere o Exemplo 5.1.  
(a) Indique os valores de $n$ e $p$.  
(b) No modelo de regressão linear simples atribuído a `fit`, indique os números das equações para o cálculo das medidas `Estimate`, `Std. Error` e `t value`.  
(c) Utilizado a linguagem R, aplique as equações indicadas no item (b) de forma a obter os valores indicados em `summary(fit)`.  
(d) No modelo de regressão pela origem atribuído a `fit0`, indique os números das equações para o cálculo das medidas `Estimate`, `Std. Error` e `t value`.  
(e) Utilizado a linguagem R, aplique as equações indicadas no item (d) de forma a obter os valores indicados em `summary(fit0)`.  
(f) Qual o modelo mais indicado segundo a sua análise? Justifique.
```


## Regressão Linear Múltipla

O modelo de *regressão linear múltipla* universal/populacional é construído, pela abordagem clássica, com todos os $N$ pares ordenados do universo, e pode ser descrito pela relação a seguir.

\begin{equation}
Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \cdots + \beta_p X_{ip} + \varepsilon_i,
(\#eq:rlm)
\end{equation}
$\varepsilon_i \sim \mathcal{N}(0,\sigma_{\varepsilon})$. Por ter dimensionalidade $p$ usualmente utiliza-se notação matricial na forma

\begin{equation}
\boldsymbol{Y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon},
(\#eq:rlm2)
\end{equation}
onde $\boldsymbol{Y} = \begin{bmatrix} Y_1 \\ Y_2 \\ \vdots \\ Y_n \end{bmatrix}$, $\boldsymbol{X} = \begin{bmatrix} 1 & X_{11} & X_{12} & \cdots & X_{1p} \\ 1 & X_{21} & X_{22} & \cdots & X_{2p}  \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & X_{n1} & X_{n2} & \cdots & X_{np} \end{bmatrix}$, $\boldsymbol{\beta} = \begin{bmatrix} \beta_1 \\ \beta_2 \\ \vdots \\ \beta_n \end{bmatrix}$ e $\boldsymbol{\varepsilon} = \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{bmatrix}$.

Para a obtenção das estimativas dos parâmetros utiliza-se a Eq. \@ref(eq:beta).
\begin{equation}
\boldsymbol{\hat{\beta}} = (\boldsymbol{X'X})^{-1} (\boldsymbol{X'X}) \boldsymbol{Y}
(\#eq:beta)
\end{equation}

```{example}
No site http://archive.ics.uci.edu/ml/datasets/Energy+efficiency está disponível uma análise de energia feita por [@tsanas2012accurate] usando 12 formas diferentes de construção simuladas no [Ecotect](https://en.wikipedia.org/wiki/Autodesk_Ecotect_Analysis). Os edifícios diferem em relação à área envidraçada, à distribuição da área envidraçada e à orientação, entre outros parâmetros. Foram simuladas várias configurações como funções das características acima mencionadas para obter 768 formas de construção. O conjunto de dados detalhado a seguir compreende 768 amostras e 8 características (`X1` a `X8`), com o objetivo de prever duas respostas reais (`Y1` e `Y2`).

`X1`: Compactação Relativa  
`X2`: Superfície  
`X3`: Área da parede  
`X4`: Área do telhado  
`X5`: Altura total  
`X6`: Orientação  
`X7`: Área de Envidraçamento  
`X8`: Distribuição da Área de Envidraçamento  
`Y1`: Carga de aquecimento  
`Y2`: Carga de resfriamento  
```
```{r}
library(readxl)
url1 <- 'http://archive.ics.uci.edu/ml/machine-learning-databases/00242/ENB2012_data.xlsx'
download.file(url1, 'temp.xlsx', mode = 'wb')
energy <- read_excel('temp.xlsx')
str(energy)   # dando uma olhada nas variáveis
pairs(energy) # matriz de dispersão
fit0 <- lm(Y1 ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8, data = energy) # modelo saturado
summary(fit0)
fit1 <- step(fit0)  # filtrando as variáveis com stepwise
summary(fit1)
```
Após o primeiro ajuste atribuído a `fit0` é possível notar que o coeficiente da variável `X4` não é possível de ser calculado devido a singularidades, i.e, impossibilidade de inversão das matrizes do modelo. Sendo assim, ao modelo saturado (contendo todas as variáveis candidatas) aplicou-se o método de *stepwise*, proposto por [@efroymson1960multiple] e utilizado para selecionar variáveis. Este método busca automaticamente o melhor conjunto de variáveis de maneira a minimizar alguma medida, usualmente o *Critério de Informação de Akaike* (AIC, na sigla em inglês), sugerido por [@akaike1974new]. De acordo com a métrica do stepwise, quanto menor o valor de AIC, melhor a combinação das variáveis.  

Pelos resultados obtidos pode-se verificar que o modelo ajustado em `fit1` possui todas as variávies significativas para um $\alpha$ inferior a 0.01, estatística F de 1387 para 6 e 761 graus de liberdade, levando a um p-value geral do modelo menor que $2.2 \times 10^{-16}$, o que indica boa aderência aos dados. O valor do `Multiple R-squared` é de 0.9162, indicando que o modelo explica em torno de 92% da variação de `Y1`. Desta forma este é um modelo aceitável, que possui coeficientes de `X1` e `X2` negativos, indicando que um aumento destas variáveis (respectivamente compactação relativa e superfície) deve reduzir a carga de aquecimento. Mais especificamente, um aumento de uma unidade na compactação relativa (`X1`) gera uma redução esperada de 64.77 unidades na carga de aquecimento, mantidas constantes as demais variáveis. As variáveis `X3`, `X5`, `X7` e `X8` possuem coeficientes positivos, levando a um impacto esperado positivo em `Y1`. Como exemplo, para cada aumento de uma unidade na altura total (`X5`) espera-se um aumento aproximando de 4.17 unidades na carga de aquecimento. As outras variáveis possuem interpretação análoga, devendo-se sempre observar o sinal dos coeficientes.

```{exercise}
Refaça o Exemplo 5.2 utilizando a variável `Y2` como resposta, ajustando o modelo saturado, filtrando as variáveis com a função `step` e interpretando os resultados.
```

<!-- Note que boa parte dos testes de hipóteses podem ser escritos como modelos de regressão. Veja o link https://lindeloev.github.io/tests-as-linear/ para maiores detalhes. -->


##	Regressão Logística

### Variáveis binárias/dicotômicas
Em problemas aplicados é comum fazer uso de variáveis aleatórias que admitam apenas dois valores, chamadas v.a. *binárias* ou *dicotômicas*. Empresas de serviços financeiros podem estar interessados em clientes adimplentes/inadimpletes, hospitais em pacientes com/sem melhora, cientistas da computação em servidores operantes/inoperantes, etc. Começamos com um exemplo numérico para ilustrar.

```{example}
Um banco está interessado na variável aleatória $Y$: cliente inadimplente. Ela assume valor 1 se o cliente é inadimplente (sucesso) e 0 caso contrário (fracasso). Note que a terminologia sucesso/fracasso indica se a variável $Y$ foi observada (1) ou não (0), ainda que 'sucesso' possa não significar algo agradável.  
Se considerarmos 40 clientes, 20 inadimplentes e 20 adimplentes, pode-se calcular a probabilidade (incondicional) de observar um ciente inadimplente (sucesso) por $Pr(Y=1)=\frac{20}{40}=0.5$. Da mesma forma, a probabilidade (incondicional) de observar um ciente adimplente (fracasso) é dada por $Pr(Y=0)=\frac{20}{40}=0.5$. Utilizando a linguagem R pode-se gerar a sequência de zeros e uns descrita acima, bem como as respectivas probabilidades.
```{r}
n <- 20
(y <- c(rep(0,n), rep(1,n)))
(taby <- table(y))
prop.table(taby)
```

Vendo a questão desta maneira pode-se considerar que a probablidade de um cliente ser inadimplente é de 50%. É possível, porém, considerar outras variáveis para refinar esta probabilidade. Suponha uma variável $X_1$, que ocorre em aproximadamente 20% dos clientes adimplentes ($Y=0$) e em aproximadamente 90% dos clientes inadimplentes ($Y=1$). Assim, visto que presença da variável $X_1$ é maior entre os inadimplentes, intuitivamente devemos atribuir uma probabilidade maior de inadimplência a clientes que apresentem a característica $X_1$, e menor para aquele onde a característica $X_1$ é ausente. Novamente podemos fazer uso da linguagem R para gerar sequências similares àquelas consideradas na teoria.

```{r}
suppressMessages(library(tidyverse))
set.seed(1); (x1 <- c(rbinom(n,1,.2), rbinom(n,1,.9)))  # gerando sequência pseudoaleatória
```

Note na saída abaixo que o segmento indicado por $\texttt{\$}$&#96;$\texttt{0}$&#96; refere-se ao grupo onde $X_1=0$, e $\texttt{\$}$&#96;$\texttt{1}$&#96; onde $X_1=1$.
```{r}
by(y,x1,table) %>%   # contando o número de zeros e uns em y separado por x1
  lapply(prop.table) # calculando as proporções 
```
Note que a probabilidade de inadimplência no grupo onde $X_1$ não é observada é igual a $Pr(Y=1|X_1=0) = 0.05882353$. No grupo de clientes que apresentam a característica $X_1$  esta probabilidade sobe para $Pr(Y=1|X_1=1) = 0.826087$.  

Finalmente é considerada uma variável $X_2$, que aparece em aproximadamente metade dos clientes inadimplentes e em aproximadamente metade dos adimplentes. Assim, observar a característica $X_2$ não deve trazer informação sobre a probabilidade de inadimplência, simbolizada pela probabilidade condicional $Pr(Y=1|X_2 = x_2)$, $x_2 \in \{0,1\}$.

```{r}
set.seed(1); (x2 <- rbinom(2*n,1,.5)) # simulando valores de X2
by(y,x2,table) %>%  # contando os zeros e uns de y, separados pelos valores de X2
  lapply(prop.table) # transformando a contagem em proporção
```

### O modelo de regressão logística
A *regressão logística* pertence à classe dos *modelos lineares generalizados*, descrita em detalhes por [@mccullagh1989generalized], [@agresti2007introduction] e [@paula2013modelos]. Seja $Y$ uma variável aleatória binária com distribuição binimial de probablidade de sucesso $\pi(x)$. A notação $\pi(x)$ sugere que a probabilidade de sucesso está condicionada a um valor/categoria $x$. Desta forma, $\pi(x) = Pr(Y=1|X=x)$. Define-se a função *logito* conforme a Eq. \@ref(eq:logit_uni), onde $\mathrm{log}$ indica o logaritmo na base $e \approx 2.718281828459$.

\begin{equation}
\mathrm{logito}\left[ \pi(x) \right] = \mathrm{log} \left( \dfrac{\pi(x)}{1-\pi(x)} \right) = \beta_0 + \beta_1 x
(\#eq:logit_uni)
\end{equation}
Isolando $\pi(x)$ na Eq. \@ref(eq:logit_uni) obtém-se

\begin{equation}
\pi(x) = \dfrac{e^{\beta_0 + \beta_1 x}}{1+e^{\beta_0 + \beta_1 x}}
(\#eq:pix_uni)
\end{equation}

```{example}
Considerando novamente as informações do Exemplo 8.1, vamos agora utilizar a estrutura das Equações \@ref(eq:logit_uni) e \@ref(eq:pix_uni) para abordar o problema.
```

```{r}
# modelo 1
x1 <- as.factor(x1)  # convertendo em fator para usar a função glm
fit1 <- glm(y ~ x1, family = 'binomial') # ajustando modelo logístico
summary(fit1) # detalhamento do modelo
```

Note que o intercepto $\hat{\beta}_0 = -2.773$ é significativo ($p-value = 0.00715$), bem como o coeficiente que indica a presença do atributo $X_1$, $\hat{\beta}_1 = 4.331$ ($p-value = 0.00021$). Assim, pode-se considerar o modelo conforme Eq. \@ref(eq:logit_uni), na forma \[\mathrm{log} \left( \dfrac{\pi(x)}{1-\pi(x)} \right) = -2.773 + 4.331 I_{x_1}. \] A simbologia $I_{x_1}$ indica uma variável *indicadora* da presença do atributo $X_1$. Assim, se a pessoa possuir o atributo $X_1$ considera-se $I_{x_1}=1$, e $I_{x_1}=0$ caso contrário. Este é um modelo conhecido como *casela de referência*, visto que a variável $X_1$ é categórica. Desta forma, uma das categorias/níveis da variável é escolhida como a (casela de) referência. Por padrão, o R utiliza a primeira da ordem numérica/alfabética, no caso $X_1=0$. Desta forma, se a pessoa não possui a característica $X_1$, pode-se calcular sua probabilidade de inadimplência pela Eq. \@ref(eq:pix_uni), dada por \[Pr(Y=1|X_1=0) = \pi(0) = \dfrac{e^{-2.773 + 4.331 \times 0}}{1+e^{-2.773 + 4.331 \times 0}} \approx 0.05882353.\] No caso de alguém que possui a característica $X_1$, \[Pr(Y=1|X_1=1) = \pi(1) = \dfrac{e^{-2.773 + 4.331 \times 1}}{1+e^{-2.773 + 4.331 \times 1}} \approx 0.826087.\]

```{r}
# modelo 1
exp(coef(fit1)[1])/(1+exp(coef(fit1)[1])) # Pr(Y = 1 | X1 = 0)
exp(sum(coef(fit1)))/(1+exp(sum(coef(fit1)))) # Pr(Y = 1 | X1 = 1)
```

O mesmo procedimento considerado para $X_1$ pode ser realizado para $X_2$. Note que o intercepto $\hat{\beta}_0 = 0.2877$ é **não** significativo ($p-value = 0.514$), bem como o coeficiente que indica a presença do atributo $X_2$, $\hat{\beta}_1 = -0.6061$ ($p-value = 0.344$).

```{r}
# modelo 2
x2 <- as.factor(x2)  # convertendo em fator para usar a função glm
fit2 <- glm(y ~ x2, family = 'binomial') # ajustando modelo logístico
summary(fit2) # detalhamento do modelo
```

Desta maneira o modelo não deve ser utilizado, mas para efeito de comparação com os resultados do Exemplo 8.1 são calculadas as probabilidades de sucesso condicionadas a $X_2=0$ e $X_2=1$.

```{r}
# modelo 2
exp(coef(fit2)[1])/(1+exp(coef(fit2)[1])) # Pr(Y = 1 | X2 = 0)
exp(sum(coef(fit2)))/(1+exp(sum(coef(fit2)))) # Pr(Y = 1 | X2 = 1)
```


```{exercise}
Considere o conjunto de dados apresentado por Ronny Kohavi e Barry Becker, disponível em https://archive.ics.uci.edu/ml/datasets/adult. Considere um modelo de regressão logística para avaliar as características que mais impactam no salario das pessoas (acima ou abaixo de 50 mil dólares).
```

```{r}
# lendo e arrumando os dados
dat <- read.table('https://archive.ics.uci.edu/ml/machine-learning-databases/adult/adult.data',
                  sep = ',')
dat <- dat[,-c(3,9)]
colnames(dat) <- c('idade','tipoTrabalho','educacao','anosEstudo','estadoCivil','ocupacao',
                   'relacao','genero','ganhoCapital','perdaCapital','horasPorSemana',
                   'paisOrigem','salario')
```