lake_huron.Rmd

---
title: "LFO-CV of Lake Huron Data"
author: "Paul Bürkner & Aki Vehtari"
output:
  html_document:
    theme: default
encoding: UTF-8
editor_options: 
  chunk_output_type: console
---

<style>
pre {
  overflow-x: auto;
}
pre code {
  word-wrap: normal;
  white-space: pre;
}
</style>

---

# {.tabset}

## Setup

```{r setup, include=FALSE, cache=FALSE}
knitr::opts_chunk$set(cache = TRUE, fig.width = 10)
```

```{r general_options, cache=FALSE}
source("sim_functions.R")
library(tidyverse)
library(brms)
options(width = 160, mc.cores = 2)
theme_set(theme_default())
```

```{r}
data("LakeHuron")
N <- length(LakeHuron)
df <- data.frame(
  y = as.numeric(LakeHuron),
  year = as.numeric(time(LakeHuron)),
  time = 1:N
)
```

## Model

Fit an AR(4) model:
```{r}
fit_lh <- brm(
  y | mi() ~ 1, 
  data = df, 
  autocor = cor_ar(~time, p = 4), 
  prior = prior(normal(0, 0.5), class = "ar"),
  chains = 2, warmup = 1000, iter = 5000,
  control = list(adapt_delta = 0.99),
  seed = 5838296, file = "models/fit_lh"
)
summary(fit_lh)
```

```{r}
preds <- posterior_predict(fit_lh)
preds <- cbind(
  Estimate = colMeans(preds), 
  Q5 = apply(preds, 2, quantile, probs = 0.05),
  Q95 = apply(preds, 2, quantile, probs = 0.95)
)

ggplot(cbind(df, preds), aes(x = year, y = Estimate)) +
  geom_smooth(aes(ymin = Q5, ymax = Q95), stat = "identity", size = 0.5) +
  geom_point(aes(y = y)) + 
  labs(y = "Water Level (ft)", x = "Year")
```

## LFO-CV of 1-SAP {.tabset}

```{r}
L <- 20
M <- 1
```

### ELPD

```{r, results="hide", warning=FALSE, message=FALSE}
lh_elpd_1sap_exact <- compute_lfo(
  fit_lh, type = "exact", M = M, L = L, 
  file = "results/lh_elpd_1sap_exact.rds"
)
lh_elpd_1sap_approx <- compute_lfo(
  fit_lh, type = "approx", M = M, L = L, 
  file = "results/lh_elpd_1sap_approx.rds"
)
```

```{r}
summarize_elpds(lh_elpd_1sap_exact)
summarize_elpds(lh_elpd_1sap_approx)
```

```{r}
ks <- na.omit(attributes(lh_elpd_1sap_approx)$ks)
ids <- (L + 1):N
plot_ks(ks, ids)
```

Perform approximate LOO-CV:
```{r}
(loo_lh <- loo(fit_lh, newdata = df[-seq_len(L), ]))
```

### RMSE

```{r, results="hide", warning=FALSE, message=FALSE}
lh_rmse_1sap_exact <- compute_lfo(
  fit_lh, type = "exact", M = M, L = L, criterion = "rmse",
  file = "results/lh_rmse_1sap_exact.rds"
)
lh_rmse_1sap_approx <- compute_lfo(
  fit_lh, type = "approx", M = M, L = L, criterion = "rmse",
  file = "results/lh_rmse_1sap_approx.rds"
)
```

```{r}
summarize_elpds(lh_rmse_1sap_exact)
summarize_elpds(lh_rmse_1sap_approx)
```


## LFO-CV of 4-SAP {.tabset}

```{r}
L <- 20
M <- 4
```

### ELPD

```{r, results="hide", warning=FALSE, message=FALSE}
lh_elpd_4sap_exact <- compute_lfo(
  fit_lh, type = "exact", M = M, L = L, 
  file = "results/lh_elpd_4sap_exact.rds"
)
lh_elpd_4sap_approx <- compute_lfo(
  fit_lh, type = "approx", M = M, L = L, 
  file = "results/lh_elpd_4sap_approx.rds"
)
```

```{r}
summarize_elpds(lh_elpd_4sap_exact)
summarize_elpds(lh_elpd_4sap_approx)
```

### RMSE

```{r, results="hide", warning=FALSE, message=FALSE}
lh_rmse_4sap_exact <- compute_lfo(
  fit_lh, type = "exact", M = M, L = L, criterion = "rmse",
  file = "results/lh_rmse_4sap_exact.rds"
)
lh_rmse_4sap_approx <- compute_lfo(
  fit_lh, type = "approx", M = M, L = L, criterion = "rmse",
  file = "results/lh_rmse_4sap_approx.rds"
)
```

```{r}
summarize_elpds(lh_rmse_4sap_exact)
summarize_elpds(lh_rmse_4sap_approx)
```

## LFO-CV of block-1-SAP {.tabset}

```{r}
L <- 20
M <- 1
B <- 10
```

### ELPD

```{r, results="hide", warning=FALSE, message=FALSE}
lh_elpd_block1sap_exact <- compute_lfo(
  fit_lh, type = "exact", M = M, L = L, B = B,
  file = "results/lh_elpd_block1sap_exact.rds"
)
lh_elpd_block1sap_approx <- compute_lfo(
  fit_lh, type = "approx", M = M, L = L, B = B,
  file = "results/lh_elpd_block1sap_approx.rds"
)
```

```{r}
summarize_elpds(lh_elpd_block1sap_exact)
summarize_elpds(lh_elpd_block1sap_approx)
```

```{r}
ks <- na.omit(attributes(lh_elpd_block1sap_approx)$ks)
ids <- (L + 1):N
plot_ks(ks, ids)
```

### RMSE

```{r, results="hide", warning=FALSE, message=FALSE}
lh_rmse_block1sap_exact <- compute_lfo(
  fit_lh, type = "exact", M = M, L = L, B = B, criterion = "rmse",
  file = "results/lh_rmse_block1sap_exact.rds"
)
lh_rmse_block1sap_approx <- compute_lfo(
  fit_lh, type = "approx", M = M, L = L, B = B, criterion = "rmse",
  file = "results/lh_rmse_block1sap_approx.rds"
)
```

```{r}
summarize_elpds(lh_rmse_block1sap_exact)
summarize_elpds(lh_rmse_block1sap_approx)
```


## LFO-CV of block-4-SAP {.tabset}

```{r}
L <- 20
M <- 4
B <- 10
```

### ELPD

```{r, results="hide", warning=FALSE, message=FALSE}
lh_elpd_block4sap_exact <- compute_lfo(
  fit_lh, type = "exact", M = M, L = L, B = B,
  file = "results/lh_elpd_block4sap_exact.rds"
)
lh_elpd_block4sap_approx <- compute_lfo(
  fit_lh, type = "approx", M = M, L = L, B = B,
  file = "results/lh_elpd_block4sap_approx.rds"
)
```

```{r}
summarize_elpds(lh_elpd_block4sap_exact)
summarize_elpds(lh_elpd_block4sap_approx)
```

### RMSE

```{r, results="hide", warning=FALSE, message=FALSE}
lh_rmse_block4sap_exact <- compute_lfo(
  fit_lh, type = "exact", M = M, L = L, B = B, criterion = "rmse",
  file = "results/lh_rmse_block4sap_exact.rds"
)
lh_rmse_block4sap_approx <- compute_lfo(
  fit_lh, type = "approx", M = M, L = L, B = B, criterion = "rmse",
  file = "results/lh_rmse_block4sap_approx.rds"
)
```

```{r}
summarize_elpds(lh_rmse_block4sap_exact)
summarize_elpds(lh_rmse_block4sap_approx)
```