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

```{r setup}
#| include: false
library(tidyverse)
library(patchwork)
library(dukestm)

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

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

```{r}
#| include: false
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
  )
```


## `dukestm` package

This is a companion package for the course where I will be putting useful functions
for some of our common tasks in the course.

This is not a official or polished package but I've tried to include documentation for most functions. If you notice a problem open and issue or send me an email.

To install,

```{r}
#| eval: false
devtools::install_github("sta344-644-fa23/dukestm")
```

<br/>

Updates will be made throughout the semester and I will attempt to remind you when something new is available, reruning the above will get you the latest version.

## Example - `epred_draws_fix()`

```{r}
#| include: false
options(width=50)
```

:::: {.columns .small}
::: {.column width='50%'}
```{r}
tidybayes::epred_draws(b, newdata=d)
```
:::

::: {.column width='50%' .fragment}
```{r}
dukestm::epred_draws_fix(b, newdata=d)
```
:::
::::

```{r}
#| include: false
options(width=80)
```

::: {.aside}
`dukestm` also has implementations of `predicted_draws_fix()` and `residual_draws_fix()`
:::


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

::: {.medium}
```{r}
( dukestm::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?

::: {.small}
```{r}
epred_draws_fix(b, newdata=d) |>
  ggplot(aes(x=x)) +
    ggdist::stat_interval(alpha=0.3, aes(y=.epred, group=x), linewidth=1.66) +
    geom_point(data=d, aes(y=y))
```
:::

## The right predictions

::: {.small}
```{r}
predicted_draws_fix(b, newdata=d) |>
  ggplot(aes(x=x)) +
    ggdist::stat_interval(alpha=0.3, aes(y=.prediction, group=x), linewidth=1.25) +
    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 -$\hat{y}$ vs $y$

:::: {.columns .small}
::: {.column width='50%'}
```{r}
y_hat_rmse = epred_draws_fix(b, newdata=d) |> 
  group_by(.iteration, .chain) |>
  yardstick::rmse(truth = y, estimate = .epred)
```

```{r}
summarize(y_hat_rmse, mean(.estimate), .by = .chain)
```

```{r}
#| echo: false
y_hat_rmse |>
  ggplot(aes(x=.estimate, fill=as.factor(.chain))) +
    geom_density(alpha=0.2) +
    guides(fill="none")
```
:::

::: {.column width='50%'}
```{r}
y_rmse = predicted_draws_fix(b, newdata=d) |> 
  group_by(.iteration, .chain) |>
  yardstick::rmse(truth = y, estimate = .prediction)
```

```{r}
summarize(y_rmse, mean(.estimate), .by = .chain)
```

```{r}
#| echo: false
y_rmse |>
  ggplot(aes(x=.estimate, fill=as.factor(.chain))) +
    geom_density(alpha=0.2) +
    guides(fill="none")
```
:::
::::



## CRPS

:::: {.columns .small}
::: {.column width='50%'}
$\hat{y}$
```{r}
epred_draws_fix(b, newdata=d) |> 
  group_by(.chain, x) |>
    summarise(
      crps = dukestm::calc_crps(.epred, obs=y), .groups="drop_last"
    ) |>
    summarize(
      mean(crps)
    )
```
:::

::: {.column width='50%'}
$y$
```{r}
predicted_draws_fix(b, newdata=d) |> 
  group_by(.chain, x) |>
    summarise(
      crps = dukestm::calc_crps(.prediction, obs=y), .groups="drop_last"
    ) |>
    summarize(
      mean(crps)
    )
```
:::
::::

## Posterior sampling functions

<br/>

![](imgs/Lec04/normal_heiss.png){fig-align="center" width="100%"}


::: {.aside}
From Andrew Heiss' blog post [Visualizing the differences between Bayesian posterior predictions, linear predictions, and the expectation of posterior predictions]( https://www.andrewheiss.com/blog/2022/09/26/guide-visualizing-types-posteriors/)
:::

# 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 (as a factor) ...

. . . 

```{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 a 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")
```

## Comparing Residuals - scale free

```{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, scales = "free_y") + 
    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")
```

## Quadratic fit - residuals (scale free)

```{r}
#| echo: false
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, scales = "free_y") + 
    guides(color="none")
```


# Bayesian Model

## Bayesian Poisson Regression Model


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

## Model priors

::: {.medium}
```{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

::: {.small}
```{r}
#| output-location: column
aids_base +
  ggdist::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

::: {.small}
```{r}
#| output-location: column
aids_base +
  ggdist::stat_lineribbon(
    data = tidybayes::predicted_draws(
      g_bayes,
      newdata = tibble(year=seq(1981,1993,by=0.1))
    ),
    aes(y=.prediction),
    alpha=0.25
  )
```
:::

## Residuals

```{r}
( g_bayes_resid = residual_draws_fix(
  g_bayes, newdata = aids
) )
```

## Residual plot

::: {.small}
```{r}
g_bayes_resid |> 
  ggplot(aes(y = .residual, x = year, color=as.factor(.chain), group=.chain)) +
    ggdist::stat_pointinterval(alpha=0.2) +
    geom_hline(yintercept = 0, color='grey', linetype=2)
```
:::

## Standardized residuals?

::: {.small}
```{r}
residual_draws_fix(
  g_bayes, newdata = aids, type = "pearson"
) |> 
  ggplot(aes(y = .residual, x = year, color=as.factor(.chain), group=.chain)) +
    ggdist::stat_pointinterval(alpha=0.2) +
    geom_hline(yintercept = 0, color='grey', linetype=2)
```
:::

## Model performance - rmse, crps

::: {.medium}
```{r}
predicted_draws_fix(g_bayes, newdata = aids) |>
  group_by(.chain, .row) |>
  summarize(
    rmse = yardstick::rmse_vec(cases, .prediction),
    crps = dukestm::calc_crps(.prediction, cases)
  ) |>
  group_by(.chain) |>
  summarize(
    rmse = mean(rmse),
    crps = mean(crps)
  )
```
:::

## Model performance - emp coverage

::: {.medium}
```{r}
predicted_draws_fix(g_bayes, newdata = aids) |>
  group_by(.row, cases) |>
  ggdist::mean_hdi(
    .prediction, .width = c(0.5,0.9,0.95)
  ) |>
  mutate(contains = cases >= .lower & cases <= .upper) %>%
  group_by(.width) |>
  summarize(
    emp_cov = sum(contains)/n()
  )
```
:::