---
title: "Diagnostics and<br/>Model Evaluation"
subtitle: "Lecture 03"
author: "Dr. Colin Rundel"
footer: "Sta 344 - Fall 2022"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
execute:
  echo: true
---

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

#source("../util-crps.R")

knitr::opts_chunk$set(
  collapse = TRUE,
  fig.width=7,
  fig.height=3.8,
  out.width="\\textwidth",
  fig.align="center",
  echo=TRUE,
  warning=FALSE
)

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


# Some more linear models

## Linear model and data

```{r echo=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
)
```

```{r}
#| fig.height: 4
ggplot(d, aes(x=x,y=y)) + 
  geom_point() + 
  geom_smooth(method="lm", color="blue", se = FALSE)
```

## Linear model

::: {.small}
```{r}
#| code-line-numbers: "|13,14,18"
l = lm(y ~ x, data=d)
summary(l)
```
:::


## Bayesian model (brms)

::: {.small}
```{r}
#| message: false
#| warning: false
#| code-line-numbers: "|20,21,25"
( 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
  )
)
```
:::

## Parameter estimates

```{r}
plot(b)
```

## tidybayes - gather_draws (long)

```{r echo=TRUE}
b_post = b |>
  tidybayes::gather_draws(b_Intercept, b_x, sigma)

b_post
```

## tidybayes - spread_draws (wide)

```{r echo=TRUE}
b_post_wide = b |>
  tidybayes::spread_draws(b_Intercept, b_x, sigma)

b_post_wide
```


## Posterior plots

```{r echo=TRUE, fig.height=4}
b_post |>
  ggplot(aes(fill=as.factor(.chain), group=.chain, x=.value)) +
  geom_density(alpha=0.33, color=NA) +
  facet_wrap(~.variable, scales = "free")
```

## Trace plots

```{r echo=TRUE, fig.height=3.8} 
b_post %>% 
  ggplot(aes(x=.iteration, y=.value, color=as.factor(.chain))) +
  geom_line(alpha=0.5) +
  facet_grid(.variable~.chain, scale="free_y") +
  geom_smooth(method="loess") + labs(color="chain")
```


## Credible Intervals

::: {.medium}
```{r echo=TRUE}
( b_ci = b_post |>
    group_by(.chain, .variable) |>
    tidybayes::mean_hdi(.value, .width=c(0.8, 0.95))
)
```
:::

## Aside - `mean_qi()` vs `mean_hdi()`

These differ in the use of the quantile interval vs. the highest-density interval.

```{r}
#| echo: false

set.seed(1234)
ci_ex = data_frame(
  dist_1 = rnorm(10000),
  dist_2 = c(rnorm(5000, 2), rnorm(5000, -2))
) %>% 
  mutate_all(function(x) (x-min(x))/max(x) ) %>%
  tidyr::gather(dist, x) 

tmp = rbind(
  ci_ex %>% group_by(dist) %>% tidybayes::mean_qi( x, .width=c(0.5,0.95)) %>% mutate(method="qi"),
  ci_ex %>% group_by(dist) %>% tidybayes::mean_hdi(x, .width=c(0.5,0.95)) %>% mutate(method="hpd")
) %>%
  group_by(method,dist, .width) %>%
  mutate(region = 1:n()) %>%
  tidyr::gather(type, bound, .lower, .upper)

tmp_range = tmp |>
  pivot_wider(c(dist, .width, method, region), names_from = type, values_from = bound)

ggplot(tmp) +
  geom_density(data=ci_ex, aes(x=x, group=dist), fill="grey", color="black") + 
  geom_rect(
    data = tmp_range,
    aes(xmin=.lower, xmax=.upper, fill=as.factor(.width), group=region), 
    ymin=-Inf, ymax=Inf, color=NA,
    alpha=0.20
  ) +
  geom_vline(aes(color=as.factor(.width), xintercept=bound), size=1.5) +
  facet_grid(method~dist, scale="free_x") +
  ylim(0,4.75) +
  labs(color="prob") +
  guides(fill="none")
```

## Caterpillar Plots

```{r echo=TRUE, fig.height=4}
b_ci %>%
  ggplot(aes(x=.value, y=.chain, color=as.factor(.chain))) + 
  facet_grid(.variable ~ .) +
  tidybayes::geom_pointintervalh() +
  ylim(0.5, 4.5)
```


# Predictions

## lm predictions

::: {.medium}
```{r}
(l_pred = broom::augment(l, interval="confidence"))
```
:::

::: {.aside}
There are the same values you would get from `predict(l, se.fit = TRUE)` just combined with `d`
:::

## Confidence interval

```{r}
l_pred |>
  ggplot(aes(x=x,y=y)) + 
   geom_point() +
   geom_line(aes(y=.fitted), col="blue", size=1.5) +
   geom_ribbon(aes(ymin=.lower, ymax=.upper), col="blue", fill="blue", alpha=0.33)
```


## brms predictions

```{r}
(b_pred = predict(b))
```

## Credible interval

```{r}
d |>
  bind_cols(b_pred) |>
  ggplot(aes(x=x,y=y)) + 
    geom_point() +
    geom_line(aes(y=Estimate), col="red", size=1.5) +
    geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), col='red', fill='red', alpha=0.33)
```


# Why are the intervals different?


## Raw predictions

```{r}
dim( brms::posterior_predict(b) )
dim( brms::posterior_epred(b) )
```


## Tidy raw predictions

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

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



## Posterior predictions vs Expected posterior predictions

::::: {.small}
:::: {.columns}
::: {.column width='50%'}
```{r echo=TRUE}
b_post_pred |>
  filter(.draw <= 25) |>
  ggplot(aes(x=x,y=y)) +
    geom_point() +
    geom_line(
      aes(y=.prediction, group=.draw), 
      alpha=0.33
    )
```
:::

::: {.column width='50%'}
```{r echo=TRUE}
b_post_epred |>
  filter(.draw <= 25) |>
  ggplot(aes(x=x,y=y)) +
    geom_point() +
    geom_line(
      aes(y=.epred, group=.draw), 
      alpha=0.33
    )
```
:::
::::
:::::


## Credible intervals

::::: {.small}
:::: {.columns}
::: {.column width='50%'}
```{r echo=TRUE}
b_post_pred |>
  ggplot(aes(x=x, y=y)) +
    geom_point() +
    tidybayes::stat_lineribbon(
      aes(y=.prediction), alpha=0.25
    )
```
:::

::: {.column width='50%'}
```{r echo=TRUE}
b_post_epred |>
  ggplot(aes(x=x, y=y)) +
    geom_point() +
    tidybayes::stat_lineribbon(
      aes(y=.epred), alpha=0.25
    )
```
:::
::::
:::::

## Posterior predictive checks

```{r}
brms::pp_check(b, ndraws = 25)
```

  
## Residuals - lm

```{r}
l |>
  broom::augment() |>
  ggplot(aes(x=x, y=.resid)) +
    geom_point() +
    geom_hline(yintercept=0, color='grey', linetype=2)
```

  
## Residual posteriors - brms

```{r}
b %>%
  tidybayes::residual_draws(newdata=d) |>
  ggplot(aes(x=x, y=.residual, group=x)) +
    tidybayes::stat_pointinterval() +
    geom_hline(yintercept = 0, color='grey', linetype=2)
```


# Model Evaluation

## Model assessment

If we remember back to our first regression class, one common option is $R^2$ which gives us the variability in $y$ explained by our model.

Quick review:
  
. . . 

$$ 
\underset{\text{Total}}{\sum_{i=1}^n \left(y_i - \bar{y}\right)^2} = \underset{\text{Model}}{\sum_{i=1}^n \left(\hat{y}_i - \bar{y}\right)^2} + \underset{\text{Error }}{\sum_{i=1}^n \left(y_i - \hat{y}_i\right)^2}
$$
  
. . .

$$
R^2 
  = \frac{SS_{model}}{SS_{total}}
  = \frac{\sum_{i=1}^n \left(\hat{Y}_i - \bar{Y}\right)^2}{\sum_{i=1}^n \left(Y_i - \bar{Y}\right)^2}
  = \frac{\text{Var}(\hat{\boldsymbol{Y}}) }{ \text{Var}({\boldsymbol{Y}}) }
  = \text{Cor}(\boldsymbol{Y}, \hat{\boldsymbol{Y}})^2 
$$

## Some data prep

::: {.medium}
```{r}
#| code-line-numbers: "|3|4-7|"
( 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
    )
)
```
:::

## Bayesian $R^2$

::: {.small}
When we compute any statistic for our model we want to do so at each iteration so that we can obtain the posterior distribution of that particular statistic (e.g. the posterior distribution of $R^2$ in this case).
:::

. . .

::: {.small}
```{r}
#| code-line-numbers: "|4|5|8-21"
( b_R2 = b_post_full %>%
    group_by(.iteration) %>%
    summarize(
      R2_classic = var(y_hat) / var(y),
      R2_bayes   = var(y_hat) / (var(y_hat) + var(resid))
    )
)
```
:::

## Uh oh ...

```{r}
#| echo: false
b_R2 %>% 
  tidyr::pivot_longer(starts_with("R2"), names_to = "method", values_to = "R2") %>%
  ggplot(aes(x=R2, fill=method)) + 
  geom_density(alpha=0.5) +
  geom_vline(xintercept=broom::glance(l)$r.squared, size=1)
```

## Sanity check

:::: {.small}
:::: {.columns}
::: {.column width='50%'}
```{r}
( l_pred = broom::augment(l) )
```
:::

::: {.column width='50%'}
```{r}
broom::glance(l)$r.squared

var(l_pred$.fitted) / var(l_pred$y)

var(l_pred$.fitted) / (var(l_pred$.fitted) + var(l_pred$.resid))
```
:::
::::
::::


## What if we collapsed first?

Here we calculate the posterior mean of $\hat{y}$ and use that to estimate $R^2$,

```{r echo=TRUE}
b_post_full %>%
  group_by(x) %>%
  summarize(
    y_hat = mean(y_hat),
    y = mean(y),
    resid = mean(y - y_hat),
    .groups = "drop"
  ) %>%
  summarize(
    R2_classic = var(y_hat) / var(y),
    R2_bayes   = var(y_hat) / (var(y_hat) + var(resid))
  )
```


## Some problems with $R^2$

::: {.medium}
Some new issues,

* $R^2$ doesn't really make sense in the Bayesian context

  * multiple possible definitions with different properties
  
  * fundamental equality doesn't hold anymore
  
  * Possible to have $R^2 > 1$
:::

. . .

::: {.medium}
Some old issues,

* $R^2$ always increases (or stays the same) when adding a predictor

* $R^2$ is highly susceptible to over fitting

* $R^2$ is sensitive to outliers

* $R^2$ depends heavily on values of $y$ (can differ for two equivalent models)
:::

# Some Other Metrics

## Root Mean Square Error

The traditional definition of rmse  is as follows

$$ \text{RMSE} = \sqrt{ \frac{1}{n} \sum_{i=1}^n \left(y_i - \hat{y_i} \right)^2 } $$
  
. . .

In the bayesian context, we have posterior samples from each parameter / prediction of interest so we can express this as

$$  \frac{1}{m} \sum_{s=1}^{m} \sqrt{  \frac{1}{n} \sum_{i=1}^n \left(y_i - {\hat{y}}_{i,s} \right)^2 } $$

where $m$ is the number of iterations and $\hat{Y}^s_i$ is the prediction for $Y_i$ at iteration $s$.


## Continuous Rank Probability Score

::: {.small}

Another approach is the continuous rank probability score which comes from the probabilistic forecasting literature, it compares the full posterior predictive distribution to the observation / truth.

$$ \text{CRPS} = \int_{-\infty}^\infty \left(F_{\hat{y}}(z) - {1}_{z \geq y}\right)^2 dz $$

where $F_{\hat{y}}$ is the CDF of $\hat{y}$ (the posterior predictive distribution for $y$) and ${1}_{z \geq Y}$ is an indicator function which equals 1 when $z \geq y$, the true/observed value of $y$.
:::

. . .

::: {.small}
Since this calculates a score for a single probabilistic prediction we can naturally extend it to multiple predictions by calculating an average CRPS
$$
  \frac{1}{n} \sum_{i=1}^n \int_{-\infty}^\infty \left(F_{\hat{y}_i}(z) - {1}_{z \geq y_i}\right)^2 dz 
$$

:::

## CDF vs Indicator
  
```{r echo=FALSE, warning=FALSE}
d_crps_ex = data_frame(
  value = rt(10000, df=2)
)

indicator = data.frame(value=seq(-12,12,len=1000)) |>
  mutate(y = as.double(value >= 0))

ggplot(d_crps_ex, aes(x=value)) +
  geom_line(data=indicator, color="black", aes(y=y), size=1, alpha=0.5) +
  stat_ecdf(size=1, alpha=0.5, color="blue") +
  xlim(-10,10)
```

## Empirical CDF vs Indicator

```{r echo=FALSE, warning=FALSE}
d_crps_ex = data_frame(
  value = rt(10, df=2)
)

indicator = data.frame(value=seq(-12,12,len=1000)) %>% mutate(y = as.double(value >= 0))

ggplot(d_crps_ex, aes(x=value)) +
  geom_line(data=indicator, color="black", aes(y=y), size=1, alpha=0.5) +
  stat_ecdf(size=1, alpha=0.5, color="blue") +
  xlim(-5,5)
```

## Accuracy vs. Precision

```{r}
#| echo: false
#| message: false
#| warning: false
source("util-crps.R")

d_crps = tibble(
  dist1 = rnorm(10000, sd=2)+0,
  dist2 = rnorm(10000, sd=1)+0,
  dist3 = rnorm(10000, sd=2)+2,
  dist4 = rnorm(10000, sd=1)+2
) %>% 
  tidyr::gather(dist)

rmses = d_crps %>% 
  group_by(dist) %>% 
  summarise(
    rmse = (value-0)^2 %>% 
      mean() %>% 
      sqrt() %>% 
      round(3)
  )

rmse_lookup = rmses$rmse %>% setNames(rmses$dist)
rmse_labeler = function(variable, value)
  paste0(value, " (rmse = ", rmse_lookup[value],")")

g_up = ggplot(d_crps, aes(value, color=dist, fill=dist)) +
  geom_density(alpha=0.1) +
  facet_grid(~dist, labeller = rmse_labeler) + 
  geom_vline(xintercept=0)

crps = d_crps %>% 
  group_by(dist) %>% 
  summarise(crps = calc_crps(value, 0))

crps_lookup = crps$crps %>% setNames(crps$dist) %>% round(3)
crps_labeler = function(variable, value)
  paste0(value, " (crps = ", crps_lookup[value],")")

indicator = data.frame(value=seq(-10,10,len=1000)) %>% mutate(y = as.double(value >= 0))

g_low = ggplot(d_crps, aes(value, color=dist)) +
  geom_line(data=indicator, color="black", aes(y=y), size=1, alpha=0.5) +
  stat_ecdf(size=1, alpha=0.5) +
  facet_grid(~dist, labeller = crps_labeler)

gridExtra::grid.arrange(
  g_up,
  g_low
)
```

## Empirical Coverage

One final method, which assesses model calibration is to examine how well credible intervals, derived from the posterior predictive distributions of the $y$s, capture the true/observed values.

```{r}
#| echo: false
set.seed(1111)
tibble(
  x = rep(rnorm(20), 1000) + rnorm(20000),
  prediction = as.character(rep(1:20, 1000))
) %>%
  ggplot(aes(y=prediction)) +
    tidybayes::stat_interval(aes(x=x)) +
    geom_vline(xintercept = 0)
```


# Back to our example 

## RMSE - $\hat{y}$

```{r}
#| message: false
broom::augment(l) |>
  yardstick::rmse(y, .fitted)

( b_rmse = b_post_full %>%
    group_by(.chain, .iteration) %>%
    summarize(
      rmse = sqrt( sum( (y - y_hat)^2 ) / n())
    )
)
```

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

:::: {.medium}
:::: {.columns}
::: {.column width='50%'}
```{r}
b_rmse |> 
  group_by(.chain) |> 
  summarize(post_mean = mean(rmse))
```
:::

::: {.column width='50%'}
```{r}
ggplot(b_rmse) +
  geom_density(
    aes(x=rmse, fill=as.factor(.chain)), 
    alpha=0.33
  ) +
  guides(fill="none")
```
:::
::::
::::







## CRPS - $\hat{y}$

:::: {.medium}
:::: {.columns}
::: {.column width='50%'}
```{r}
#| message: false
( b_crps = b_post_full |>
    group_by(.chain, x) |>
    summarise(
      crps = calc_crps(y_hat, y)
    )
)
```
:::

::: {.column width='50%'}
```{r}
b_crps |>
  group_by(.chain) |>
  summarize( 
    avg_crps = mean(crps)
  )
```
:::
::::
::::

## Empirical Coverage - $\hat{y}$

::: {.medium}
```{r}
( b_cover = b_post_full %>% 
    group_by(x, y) %>% 
    tidybayes::mean_hdi(
      y_hat, .prob = c(0.5,0.9,0.95)
    ) 
)
```
:::

## Empirical Coverage - $\hat{y}$ - Results

::: {.medium}
```{r}
b_cover %>%
  mutate(contains = y >= .lower & y <= .upper) %>%
  group_by(prob = .width) %>%
  summarize(
    emp_cov = sum(contains)/n()
  )
```
:::

# What went wrong now?
