---
title: "Residual Analysis + Generalized Linear Models"
subtitle: "Lecture 04"
author: "Dr. Colin Rundel"
footer: "Sta 344 - Fall 2022"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
execute:
  echo: true
  warning: false
  collapse: true
---

```{r setup, include=FALSE}
library(tidyverse)
library(patchwork)

knitr::opts_chunk$set(
  fig.align = "center"
)

ggplot2::theme_set(ggplot2::theme_bw())
```


# Lecture 3 wrap up

```{r}
#| include: false
source("fix_pred_draws.R")

set.seed(01232018)
n = 100

d = tibble(
  x = 1:n,
  y = arima.sim(n=100, list(ar=0.9,sq=1)) %>% as.numeric() + x * 0.07
)

b = brms::brm(
  y ~ x, data=d,
  prior = c(
    brms::prior(normal(0, 100), class = Intercept),
    brms::prior(normal(0, 10),  class = b),
    brms::prior(cauchy(0, 2),   class = sigma)
  ),
  silent = 2, refresh = 0
)


b_post_full = b |>
  tidybayes::spread_draws(b_Intercept, b_x, sigma) |>
  tidyr::expand_grid(d) |>
  mutate(
    y_hat = b_Intercept + b_x * x,
    resid = y - y_hat
  )
```

## Where we left it - Empirical Coverage ($\hat{y}$)

::: {.medium}
```{r}
( epred_draws_fix(b, newdata=d) %>%
    group_by(x, y) %>% 
    tidybayes::mean_hdi(
      .epred, .width = c(0.5,0.9,0.95)
    ) %>%
    mutate(contains = y >= .lower & y <= .upper) %>%
    group_by(prob = .width) %>%
    summarize(
      emp_cov = sum(contains)/n()
    )
)
```
:::

## What went wrong?

```{r}
epred_draws_fix(b, newdata=d) %>%
  ggplot(aes(x=x)) +
    tidybayes::stat_interval(alpha=0.3, aes(y=.epred, group=x)) +
    geom_point(data=d, aes(y=y))
```

## The right predictions

```{r}
predicted_draws_fix(b, newdata=d) %>%
  ggplot(aes(x=x)) +
    tidybayes::stat_interval(alpha=0.3, aes(y=.prediction, group=x)) +
    geom_point(data=d, aes(y=y))
```

## Empirical Coverage ($y$)

::: {.medium}
```{r}
predicted_draws_fix(b, newdata=d) %>%
  group_by(x, y) %>% 
  tidybayes::mean_hdi(
    .prediction, .width = c(0.5,0.8,0.9,0.95)
  ) %>%
  mutate(contains = y >= .lower & y <= .upper) %>%
  group_by(prob = .width) %>%
  summarize(
    emp_cov = sum(contains)/n()
  )

```
:::

## RMSE - $y$ vs $\hat{y}$

:::: {.small}
:::: {.columns}
::: {.column width='50%'}
```{r}
epred_draws_fix(b, newdata=d) %>% 
  group_by(.iteration, .chain) %>%
  yardstick::rmse(truth = y, estimate = .epred) %>%
  ggplot(aes(x=.estimate, fill=as.factor(.chain))) +
    geom_density(alpha=0.2) +
    guides(fill="none")
```
:::

::: {.column width='50%'}
```{r}
predicted_draws_fix(b, newdata=d) %>%
  group_by(.iteration, .chain) %>%
  yardstick::rmse(truth = y, estimate = .prediction) %>%
  ggplot(aes(x=.estimate, fill=as.factor(.chain))) +
    geom_density(alpha=0.2)+
    guides(fill="none")
```
:::
::::
::::



# Residual Analysis

## Atmospheric $\text{CO}_2$ (ppm) from Mauna Loa

```{r echo=FALSE}
mauna_loa = tibble(
  co2 = c(co2), 
  date = c(time(.env$co2))
) %>% 
  mutate(
    year = floor(date),
    month = month.abb[(date %% 1)*12 + 1]
  ) %>%
  filter(year >= 1985)

co2_base = ggplot(mauna_loa, aes(x=date, y=co2)) +
  geom_line()

co2_base
```


## Where to start?

Well, it looks like stuff is going up on average ...

. . .

```{r}
l = lm(co2~date, data=mauna_loa)
```

```{r echo=FALSE}
mauna_loa_l = broom::augment(l, data=mauna_loa)

l_model = co2_base + 
  geom_line(data=mauna_loa_l, aes(y=.fitted), col='red', alpha=0.5)
l_resid = ggplot(mauna_loa_l, aes(x=date, y=.resid)) + 
  geom_point() + 
  geom_line(color="grey",size=0.5,alpha=0.5)

l_model / l_resid
```


## and then?

Well there is some periodicity lets add the month ...

. . . 

```{r}
ls = lm(.resid~month, data=mauna_loa_l)
```

```{r echo=FALSE}
mauna_loa_ls = broom::augment(ls, data=mauna_loa_l)

ls_model = l_resid + geom_line(data=mauna_loa_ls, aes(y=.fitted), col='red', alpha=0.5)
ls_resid = ggplot(mauna_loa_ls, aes(x=date, y=.resid)) + 
  geom_point() + 
  geom_line(color="grey",size=0.5,alpha=0.5)

ls_model / ls_resid
```

## and then and then?

There is still some long term trend in the data, maybe a fancy polynomial can help ...

. . . 

```{r}
lsy = lm(.resid~poly(date,5), data=mauna_loa_ls)
```

```{r echo=FALSE}
mauna_loa_lsy = broom::augment(lsy, data=mauna_loa_ls)

lsy_model = ls_resid + geom_line(data=mauna_loa_lsy, aes(y=.fitted), col='red', alpha=0.5)
lsy_resid = ggplot(mauna_loa_lsy, aes(x=date, y=.resid)) + 
  geom_point() + 
  geom_line(color="grey",size=0.5,alpha=0.5)

lsy_model / lsy_resid
```

## Putting it all together ...

::: {.small}
```{r}
l_comb = lm(co2~date + month + poly(date,5), data=mauna_loa)
summary(l_comb)
```
:::

## Combined fit + Residuals

```{r echo=FALSE, message=FALSE, warning=FALSE}
mauna_loa_comb = broom::augment(l_comb, data=mauna_loa)

lsy_model = co2_base + 
  geom_line(data=mauna_loa_comb, aes(y=.fitted), col='red', alpha=0.5)
lsy_resid = ggplot(mauna_loa_comb, aes(x=date, y=.resid)) + 
  geom_point() + 
  geom_line(color="grey",size=0.5,alpha=0.5)

lsy_model / lsy_resid
```

## Model performance

```{r}
#| echo: false
mauna_loa_lin = lm(co2~date, data=mauna_loa) %>%
  broom::augment(data=mauna_loa)
mauna_loa_month = lm(co2~month, data=mauna_loa) %>%
  broom::augment(data=mauna_loa)
mauna_loa_lin_month = lm(co2~date+month, data=mauna_loa) %>%
  broom::augment(data=mauna_loa)
mauna_loa_poly = lm(co2~poly(date,5), data=mauna_loa) %>%
  broom::augment(data=mauna_loa)
mauna_loa_poly_month = lm(co2~month+poly(date,5), data=mauna_loa) %>%
  broom::augment(data=mauna_loa)

calc_rmse = function(d) {
  yardstick::rmse(d, co2, .fitted)$.estimate %>%
    round(3)
}
```

| Model                           | rmse                                |
|:--------------------------------|:-----------------------------------:|
| `co2 ~ date`                    | `r calc_rmse(mauna_loa_lin)`        |
| `co2 ~ month`                   | `r calc_rmse(mauna_loa_month)`      |  
| `co2 ~ date+month`              | `r calc_rmse(mauna_loa_lin_month)`  |      
| `co2 ~ poly(date,5)`            | `r calc_rmse(mauna_loa_poly)`       | 
| `co2 ~ month+poly(date,5)`      | `r calc_rmse(mauna_loa_poly_month)` |       
| `co2 ~ date+month+poly(date,5)` | `r calc_rmse(mauna_loa_comb)`       | 



# Generalized Linear Models


## Background

A generalized linear model has three key components:
  
1. a probability distribution (from the exponential family) that describes your response variable

2. a linear predictor $\boldsymbol{\eta} = \boldsymbol{X}\boldsymbol{\beta}$,

3. and a link function $g$ such that $g(E(\boldsymbol{Y}|\boldsymbol{X})) = \boldsymbol{\eta}$ (or $E(\boldsymbol{Y}|\boldsymbol{X}) = g^{-1}(\boldsymbol{\eta})$).


## Poisson Regression

This is a special case of a generalized linear model for count data where we assume the outcome variable follows a poisson distribution (mean = variance).

$$
\begin{aligned}
Y_i &\sim \text{Poisson}(\lambda_i)\\
\log E(Y_i|\boldsymbol{X}_{i\cdot}) &= \log{\lambda_i} = \underset{1 \times p}{\boldsymbol{X}_{i\cdot}}\underset{p \times 1}{\boldsymbol{\beta}}
\end{aligned}
$$


## Example - AIDS in Belgium

These data represent the total number of new AIDS cases reported in Belgium during the early stages of the epidemic.

```{r echo=FALSE}
aids = tibble(
  year = 1981:1993,
  cases = c(12, 14, 33, 50, 67, 74, 123, 141, 165, 204, 253, 246, 240) %>% as.integer()
)

aids_base = ggplot(aids, aes(x=year, y=cases)) + 
  geom_point() +
  labs(title="AIDS cases in Belgium")
```

:::: {.columns}
::: {.column width='20%'}
::: {.small}
```{r}
aids
```
:::
:::
::: {.column width='80%'}
```{r}
#| echo: false
aids_base
```
:::
::::


## Frequentist glm fit

```{r}
( g = glm(cases~year, data=aids, family=poisson) )

```

## Model Fit

```{r}
g_pred = broom::augment(
  g, type.predict = "response", 
  newdata = tibble(year=seq(1981,1993,by=0.1))
)

aids_base + 
  geom_line(data=g_pred, aes(y=.fitted), size=1.2, alpha=0.3)
```

## Residuals?

The naive approach is to use standard residuals,

$$ r_i = Y_i - E(Y_i|X) = Y_i - \hat\lambda_i$$

. . .

::: {.small}
```{r}
g_pred_std = broom::augment(
  g, type.predict = "response"
) %>%
  mutate(.resid = cases - .fitted)
```
:::

```{r}
#| echo: false
#| fig-height: 3
ggplot(g_pred_std, aes(x=year, y=.resid)) + 
  geom_point() + 
  geom_segment(aes(xend=year, yend=0)) 
```


## Accounting for variability

Pearson residuals:
$$ r_i = \frac{Y_i - E(Y_i|X)}{\sqrt{Var(Y_i|X)}} = \frac{Y_i - \hat\lambda_i}{\sqrt{\hat\lambda_i}}$$

. . .


::: {.small}
```{r}
g_pred_pearson = broom::augment(
  g, type.predict = "response", type.residuals = "pearson"
)
```
:::

```{r}
#| echo: false
#| fig-height: 3
ggplot(g_pred_pearson, aes(x=year, y=.resid)) + 
  geom_point() + 
  geom_segment(aes(xend=year, yend=0)) 
```

## Deviance

Deviance is a way of measuring the difference between a GLM's fit and the fit of the perfect model (i.e. where $\theta_best = E(Y_i|X) = Y_i$).

It is defined as twice the log of the ratio between the likelihood of the perfect model and the likelihood of the given model,
$$ 
\begin{aligned}
D &= 2\log\left(
\frac{\mathcal{L}(\theta_{best}|Y)}
     { \mathcal{L}(\hat\theta|Y)}\right) \\
  &= 2\big(\mathcal{l}(\theta_{best}|Y) - \mathcal{l}(\hat\theta|Y)\big)
\end{aligned}
$$

## Derivation - Normal



## Derivation - Poisson


## glm output

```{r}
summary(g)
```


## Deviance residuals

::: {.small}
We can therefore think of deviance as $D = \sum_{i=1}^n d_i^2$ where $d_i$ is a generalized residual.

In the Poisson case we have,
$$ d_i = \text{sign}(y_i - \lambda_i) \sqrt{2(y_i \log (y_i/\hat\lambda_i) - (y_i-\hat\lambda_i))}$$
:::

. . .

::: {.small}
```{r}
g_pred_dev = broom::augment(
  g, type.predict = "response", type.residuals = "deviance"
)
```
:::


```{r}
#| echo: false
#| fig-height: 3
ggplot(g_pred_dev, aes(x=year, y=.resid)) + 
  geom_point() + 
  geom_segment(aes(xend=year, yend=0))  
```  


## Comparing Residuals

```{r}
#| echo: false
bind_rows(
  g_pred_std     %>% mutate(type = "Standard"),
  g_pred_pearson %>% mutate(type = "Pearson"),
  g_pred_dev     %>% mutate(type = "Deviance")
) %>%
  select(year, .fitted, .resid, type) %>%
  mutate(type = forcats::as_factor(type)) %>%
  ggplot(aes(x=year, y=.resid, color=type)) +
    geom_point() + geom_segment(aes(xend=year, yend=0)) +
    facet_wrap(~type) + 
    guides(color="none")
```

# Updating the model

## Quadratic fit

```{r}
g2 = glm(cases~year+I(year^2), data=aids, family=poisson)

g2_pred = broom::augment(
  g2, type.predict = "response",
  newdata=tibble(year=seq(1981,1993,by=0.1))
) 
```


```{r echo=FALSE}
aids_base + 
  geom_line(data=g2_pred, aes(y=.fitted), size=1.2, alpha=0.3)
```

## Quadratic fit - residuals

```{r}
#| echo: false
g2_pred = aids %>%
  mutate(
    .fitted = predict(g2, type="response"),
    Standard = cases - .fitted,
    Pearson  = residuals(g2, type="pearson"),
    Deviance = residuals(g2, type="deviance")
  )
  
g2_pred %>%
  tidyr::pivot_longer(c(Standard, Pearson, Deviance), names_to = "type", values_to = ".resid") %>%
  mutate(type = forcats::as_factor(type)) %>%
  ggplot(aes(x=year, y=.resid, color=type)) +
    geom_point() + geom_segment(aes(xend=year, yend=0)) +
    facet_wrap(~type) + 
    guides(color="none")
```


# Bayesian Model

## Bayesian Poisson Regression Model

```{r}
#| include: false
( g_bayes = brms::brm(
    cases~year, data=aids, family=poisson,
    silent=2, refresh=0
) )
```

::: {.small}
```{r}
#| eval: false
( g_bayes = brms::brm(
    cases~year, data=aids, family=poisson,
    silent=2, refresh=0
) )
```

```{r}
#| echo: false
g_bayes
```
:::

## Model priors

```{r}
brms::prior_summary(g_bayes)
```

## MCMC Diagnostics

```{r}
plot(g_bayes)
```

## Posterior Predictive Check

```{r}
brms::pp_check(g_bayes)
```

## Model fit - $\lambda$ CI

```{r}
aids_base +
  tidybayes::stat_lineribbon(
    data = tidybayes::epred_draws(
      g_bayes,
      newdata = tibble(year=seq(1981,1993,by=0.1))
    ),
    aes(y=.epred),
    alpha=0.25
  )
```


## Model fit - $Y$ CI

```{r}
aids_base +
  tidybayes::stat_lineribbon(
    data = tidybayes::predicted_draws(
      g_bayes,
      newdata = tibble(year=seq(1981,1993,by=0.1))
    ),
    aes(y=.prediction),
    alpha=0.25
  )
```
