---
title: "Modelos Lineares Generalizados"
author: "PET Estatística UFPR"
vignette: >  
  %\VignetteIndexEntry{Modelos Lineares Generalizados}  
  %\VignetteEngine{knitr::rmarkdown}  
  %\VignetteEncoding{UTF-8}---
---  

```{r setup, include=FALSE}
# source("config/_setup.R")
```

# Análise Exploratória dos Dados

Para exemplificar o uso de um modelo linear generalizado, também chamado de
MLG (ou GLM), vamos utilizar o conjunto de dados `PaulaEx3.7.25.txt`.   

```{r}
library(labestData)
# lista os dados dentro do pacote
# ls("package:labestData")
```
Podemos olhar para a descrição dos dados, obtida através da função `help`:
```{r, eval=FALSE}
help(PaulaEx3.7.25, help_type = "html")
```

Este banco de dados trata sobre a a pulsação em repouso de homens que 
fumam. A variável resposta é dicotômica, podendo ser pulsção "normal"
ou "alta". As duas outras variáveis, ou seja, o peso do indivíduo e o 
fator "fuma" (sim ou não), serão utilizadas como possíveis explicativas 
no nível de pulsação. 

Para cada variável, precisamos estar atentos às respectivas unidades 
de medida, pois as interpretações do ajuste do modelo dependem, também, 
disso. Aqui não teremos problemas, porque temos conhecimento e clareza 
sobre estas unidades. 

```{r}
#-----------------------------------------------------------------------
# Lendo a partir do repositório do labestData.
# url <- paste0("https://gitlab.c3sl.ufpr.br/pet-estatistica",
#               "/labestData/raw/devel/data-raw/PimentelEg5.2.txt")
## PimentelEg5.2 <-
#     read.table(file = url, sep = "\t", header = TRUE)
#-----------------------------------------------------------------------

# Iniciando a análise exploratória

# Observando a estrutura dos dados
str(PaulaEx3.7.25)

# Tabelas de frequência
xtabs(~fuma, data = PaulaEx3.7.25)
xtabs(~pulsa, data = PaulaEx3.7.25)

# Diagrama de dispersão.
# carregamento do pacote lattice
library(lattice)
xyplot(pulsa ~ peso, groups = fuma, 
    data = PaulaEx3.7.25,  type = "p",  
    auto.key = list(corner = c(0.05, 0.95),
    columns = 1, title = "Fuma", cex.title = 1.1), 
    xlab = "Peso do Indivíduo, em kg",       
    ylab = "Níveis de Pulsação em Repouso", 
    scales = list(x = list(rot = 90))
    )
    
# BoxPlots    
library(ggplot2)

ggplot(PaulaEx3.7.25, aes(x = pulsa, y = peso, fill = pulsa)) + geom_boxplot() +
    facet_wrap(~ fuma, ncol = 5) 
```

Notamos, a partir dos gráficos acima, que a maior parte dos indivíduos
que fumam têm pulsação em repouso alta. Mesmo assimEssas inferências podem
ajudar a interpretar e evidenciar os resultados do modelo. 

# Ajustando um modelo linear generalizado aos dados

Um modelo linear generalizado depende da existência de 3 componentes:

  * Um componente aleatório;
  * Um componente sistemático; e
  * Uma função de ligação, que lineariza a relação entre os componentes
    aleatório e sistemático. 
    
    
- O componente aleatório é um conjunto de variáveis aleatórias independentes
($Y_{1}, Y_{2}...Y_{n}$), com uma mesma distribuição da família exponencial de 
dispersão, cujo formato possa ser escrito da seguinte maneira:

$$ f(y_{i}; \theta_{i}, \phi) = exp \{ \frac{\theta_{i}y_{i} - b(\theta_{i})}{\phi} - c(y_{i}; \phi) \}$$

Com:

$E(Y_{i}) = b'(\theta_{i})$ e $Var(Y_{i}) = \phi b''(\theta_{i}) = \phi V(\mu_{i})$


- O componente sistemático é o preditor linear, contemplando um conjunto de 
covariáveis ($x_{1}, x_{2}...x_{p}$), por meio de uma combinação linear nos
parâmetros:
    
    $$\eta_{i} = \beta_{0} + \beta_{1}x_{i1} + \beta_{2}x_{i2} + ... + \beta_{p}x_{ip} $$
    
- A função de ligação é responsável por associar, de forma **linear** o 
componente aleatório ao sistemático. Ela é real, monótona e diferenciável. 
Sua notação é a seguinte:
    
    $$ g(\mu_{i}) = \nu_{i} = \beta_{0} + \beta_{1}x_{i1} + \beta_{2}x_{i2} + ... + \beta_{p}x_{ip}, $$
    
    e, equivalentemente, 
    
    $$ \mu_{i} = g^{-1}(\eta_{i}) = g(\beta_{0} + \beta_{1}x_{i1} + \beta_{2}x_{i2} + ... + \beta_{p}x_{ip})$$
    
Caso algum destes elementos não seja satisfeito, o modelo não é linear
generalizado.


Códigos para o ajuste:
```{r}
# como a variável 'pulsa' é um fator binário, vamos
# utilizar a distribuição binomial
m0 <- glm(pulsa ~ peso + fuma, 
          data = PaulaEx3.7.25,
          family= binomial(link="logit")) # especificação da função 
                                          # de ligação 'logito'
summary(m0) # resumo do ajuste   
levels(PaulaEx3.7.25$pulsa) # níveis da variável resposta
``` 

# Interpretando a saída do modelo

Dados os resultados acima, concluimos que a variável `peso`
contribui para uma menor incidência de pulsação alta. Isto também
ocorre para quando o indivíduo é não fumante, isto é, ele tem 
mais chances de ter uma pulsação normal em repouso. Obtemos
essas conclusões uma vez que ambos os parâmetros $\beta_{1}$ e
$\beta_{2}$ têm sinal negativo. 

Como esse é um caso *binomial*, estamo modelando uma probabilidade,
ou seja, valores que devem pertencer ao intervalo $[0,1]$. As estimativas
acima estão na escala do preditos. Para colocar na escala da **resposta**,
a seguinte equação é utilizada:

$$ \hat\pi_{x} = \frac{\exp^{\beta_{0} + \beta_{1}  x_{1} + \beta_{2} x_{2}}}{1 + \exp^{\beta_{0} + \beta_{1} x_{1} + \beta_{2} x_{2}}}$$
$$ \hat\pi_{x} = \frac{\exp^{3.18021  -0.05050  x_{1} -1.19298 x_{2}}}{1 + \exp^{3.18021  -0.05050 x_{1} -1.19298 x_{2}}}$$ 

Todavia, o valor da *deviance* residual parece ser muito alto. Dela
para a do modelo nulo, temos uma diferença de apenas 7.57. Isso 
é indicativo de má qualidade do ajuste. Logo, ao ajustar um modelo
linear generalizado, devemos atentar para este ponto, porque medidas
de melhoria do ajuste podem ser tomadas.

A análise de resíduos para um modelo linear generalizado é um pouco 
diferente da de um modelo de regressão. Neste caso, os resíduos 
dificilmente seguem uma distribuição Normal. Por isso, avaliamos
um gráfico com envelope simulado via *bootstrap*, que permite checar a
qualidade do modelo, mesmo que os resíduos não se distribuam de acordo
com uma Normal. 

Neste gráfico, espera-se que, se o modelo estiver bem ajustado, os
resíduos estarão dispersos em ordem aleatória dentro do envelope. 
Pequenas fugas devem ser analisadas levando em consideração o contexto
do modelo. Na Figura 1, temos envelope simulado para o ajuste feito 
anteriormente.

```{r, fig.cap='Análise de resíduos de um modelo linear generalizado'}
# Envelope simulado para a análise de resíduos
# carregamento do pacote 
library(hnp)
par(mfrow = c(1,1))
plot(hnp(m0, conf = 0.95), 
     main = "Envelope Simulado",
     xlab = "Quantis Teóricos",
     ylab = "Resíduos",
     pch = 16,
     col = c("black", "tomato", "black", "blue"))


```

