---
title: "Diagnostics and<br/>Model Evaluation"
subtitle: "Lecture 03"
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
  fig-align: center
---

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

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)

```{r}
#| message: false
#| warning: false
#| code-line-numbers: "|20,21,25"
#| output-location: slide
( 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}
(b_post = b |>
  tidybayes::gather_draws(b_Intercept, b_x, sigma))
```

## tidybayes - spread_draws (wide)

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


## 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}
( b_ci = b_post |>
    group_by(.chain, .variable) |>
    ggdist::mean_hdi(.value, .width=c(0.95))
)
```
:::

## `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}
#| fig-height: 4
b_post |>
  ggplot(aes(x=.value, y=.chain, color=as.factor(.chain))) + 
  facet_grid(.variable ~ .) +
  ggdist::stat_pointinterval() +
  ylim(0.5, 4.5)
```


# Predictions

## lm predictions

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

::: {.aside}
These are values you would get from `predict(l, se.fit = TRUE)` added to the data frame `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)) |> as_tibble()
```

## 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

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

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

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


## Posterior predictions vs<br/>Expected posterior predictions

::::: {.small}
:::: {.columns}
::: {.column width='50%'}
```{r}
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}
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 interval - predicted

::: {.small}
```{r}
  ggplot(b_post_pred, aes(x=x, y=y)) +
    geom_point() +
    ggdist::stat_lineribbon(
      aes(y=.prediction), alpha=0.25
    )
```
:::

## Credible interval - epred

::: {.small}
```{r}
  ggplot(b_post_epred, aes(x=x, y=y)) +
    geom_point() +
    ggdist::stat_lineribbon(
      aes(y=.epred), alpha=0.25
    )
```
:::


## Confidence interval (frequentist)

```{r}
#| echo: false
b_post_epred |>
  ggplot(aes(x=x, y=y)) +
    geom_point() +
    geom_ribbon(
      data = broom::augment(l, interval="confidence"),
      aes(ymin=.lower, ymax=.upper), col="blue", fill="blue", alpha=0.30
    ) +
    ggdist::stat_lineribbon(
      aes(y=.epred), alpha=0.25
    )
```

::: {.aside}
Frequentist interval generated by `broom::augment(l, interval="confidence")`
:::


## Prediction interval (frequentist)

```{r}
#| echo: false
b_post_pred |>
  ggplot(aes(x=x, y=y)) +
    geom_point() +
    geom_ribbon(
      data = broom::augment(l, interval="prediction"),
      aes(ymin=.lower, ymax=.upper), col="blue", fill="blue", alpha=0.30
    ) + 
    ggdist::stat_lineribbon(
      aes(y=.prediction), alpha=0.25
    )
```

::: {.aside}
Frequentist interval generated by `broom::augment(l, interval="prediction")`
:::

# Other useful plots

## 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)) +
    ggdist::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-6"
#| output-location: column
( 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)
```

## Frequentist sanity check

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

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

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


## What if we collapsed first?

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

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


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 = dukestm::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

:::: {.columns .small}
::: {.column width='50%'}
```{r}
#| message: false
broom::augment(l) |>
  yardstick::rmse(y, .fitted)
```
:::

::: {.column width='50%'}
```{r}
( b_rmse = b_post_full |>
    group_by(.chain, .iteration) |>
    summarize(
      rmse = sqrt( sum( (y - y_hat)^2 ) / n())
    )
)
```
:::
::::


## RMSE (cont)

:::: {.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

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

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

## Empirical Coverage

::: {.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?