---
title: "3. Priors"
output:
  bookdown::html_document2:
    base_format: rmarkdown::html_vignette
    toc: true
    toc_depth: 2
    number_sections: true
pkgdown:
  as_is: true
vignette: >
  %\VignetteIndexEntry{3. Priors}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  echo = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 1.5
)
```

# Set-up

```{r setup}
suppressPackageStartupMessages({
  library(bage)
  library(poputils)
  library(ggplot2)
  library(dplyr)
})	
```

```{r}
plot_exch <- function(draws) {
  ggplot(draws, aes(x = element, y = value)) +
    facet_wrap(vars(draw), nrow = 1) +
    geom_hline(yintercept = 0, col = "grey") +
    geom_point(col = "darkblue", size = 0.7) +
    scale_x_continuous(n.breaks = max(draws$element)) +
    xlab("Unit") +
    ylab("") +
    theme(text = element_text(size = 8))
}	       
```

```{r}
plot_cor_one <- function(draws) {
  ggplot(draws, aes(x = along, y = value)) +
    facet_wrap(vars(draw), nrow = 1) +
    geom_hline(yintercept = 0, col = "grey") +
    geom_line(col = "darkblue") +
    xlab("Unit") +
    ylab("") +
    theme(text = element_text(size = 8))
}	       
```


```{r}
plot_cor_many <- function(draws) {
  ggplot(draws, aes(x = along, y = value)) +
    facet_grid(vars(by), vars(draw)) +
    geom_hline(yintercept = 0, col = "grey") +
    geom_line(col = "darkblue") +
    xlab("Unit") +
    ylab("") +
    theme(text = element_text(size = 8))
}	       
```

```{r}
plot_svd_one <- function(draws) {
  ggplot(draws, aes(x = age_mid(age), y = value, color = sexgender)) +
    facet_wrap(vars(draw), nrow = 1) +
    geom_line() +
    scale_color_manual(values = c("darkgreen", "darkorange")) +
    xlab("Age") +
    ylab("") +
    theme(text = element_text(size = 8),
          legend.position = "top",
          legend.title = element_blank())
}
```


```{r}
plot_svd_many <- function(draws) {
  draws |>
    mutate(element = paste("Unit", element)) |>
    ggplot(aes(x = age_mid(age), y = value, color = sexgender)) +
    facet_grid(vars(element), vars(draw)) +
    geom_line() +
    scale_color_manual(values = c("darkgreen", "darkorange")) +
    xlab("Age") +
    ylab("") +
    theme(text = element_text(size = 8),
          legend.position = "top",
          legend.title = element_blank())
}
```

# Exchangeable Units

## Fixed Normal `NFix()`

### Model

\begin{equation}
  \beta_j  \sim \text{N}(0, \mathtt{sd}^2)
\end{equation}

### All defaults

`sd = 1`

```{r}
set.seed(0)

NFix() |>
  generate(n_element = 10, n_draw = 8) |>
  plot_exch()
```

### Reduce `sd`

`sd = 0.01`

```{r, fig.height = 1.5}
set.seed(0)

NFix(sd = 0.01) |>
  generate(n_element = 10, n_draw = 8) |>
  plot_exch()
```


## Normal `N()`

### Model

\begin{align}
  \beta_j & \sim \text{N}(0, \tau^2) \\
  \tau & \sim \text{N}^+(0, \mathtt{s})
\end{align}

### All defaults

`s = 1`

```{r}
set.seed(0)

N() |>
  generate(n_element = 10, n_draw = 8) |>
  plot_exch()
```

### Reduce `s`

`s = 0.01`

```{r, fig.height = 1.5}
set.seed(0)

N(s = 0.01) |>
  generate(n_element = 10, n_draw = 8) |>
  plot_exch()
```



# Units Correlated With Neighbours

## Random Walk `RW()`

### Model

\begin{align}
  \beta_1 & \sim \text{N}(0, \mathtt{sd}^2) \\
  \beta_j & \sim \text{N}(\beta_{j-1}, \tau^2), \quad j = 2, \cdots, J \\
  \tau & \sim \text{N}^+(0, \mathtt{s}^2)
\end{align}

### All defaults

`s = 1`, `sd = 1`

```{r}
set.seed(0)

RW() |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Reduce `s`

`s = 0.01`, `sd = 1`

```{r}
set.seed(0)

RW(s = 0.01) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Reduce `s` and `sd`


`s = 0.01`, `sd = 0`

```{r}
set.seed(0)

RW(s = 0.01, sd = 0) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```


## Second-Order Random Walk `RW2()`

### Model

\begin{align}
  \beta_1 & \sim \text{N}(0, \mathtt{sd}^2) \\
  \beta_2 & \sim \text{N}(\beta_1, \mathtt{sd\_slope}^2) \\
  \beta_j & \sim \text{N}(2\beta_{j-1} - \beta_{j-2}, \tau^2), \quad j = 3, \cdots, J \\
  \tau & \sim \text{N}^+(0, \mathtt{s}^2)
\end{align}

### All defaults

`s = 1`, `sd = 1`, `sd_slope = 1`

```{r}
set.seed(0)

RW2() |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Reduce `s`

`s = 0.01`, `sd = 1`, `sd_slope = 1`

```{r}
set.seed(0)

RW2(s = 0.01) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```


### Reduce `s` and `sd_slope`

`s = 0.01`, `sd = 1`, `sd_slope = 0.01`

```{r}
set.seed(0)

RW2(s = 0.01, sd_slope = 0.01) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Reduce `s`, `sd`, and `sd_slope`

`s = 0.01`, `sd = 0`, `sd_slope = 0.01`

```{r}
set.seed(0)

RW2(s = 0.01, sd = 0, sd_slope = 0.01) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```


## Autoregressive `AR()`

### Model

\begin{equation}
  \beta_j \sim \text{N}\left(\phi_1 \beta_{j-1} + \cdots + \phi_{\mathtt{n\_coef}} \beta_{j-\mathtt{n\_coef}},  \omega^2\right)
\end{equation}
TMB derives a value of $\omega$ that gives each $\beta_j$ variance $\tau^2$. The prior for $\tau$ is
\begin{equation}
  \tau  \sim \text{N}^+(0, \mathtt{s}^2).
\end{equation}
The prior for each $\phi_k$ is
\begin{equation}
  \frac{\phi_k + 1}{2} \sim \text{Beta}(\mathtt{shape1}, \mathtt{shape2}).
\end{equation}

### All defaults

`n_coef = 2`, `s = 1`

```{r}
set.seed(0)

AR() |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Increase `n_coef`

```{r}
set.seed(0)

AR(n_coef = 3) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```


### Reduce `s`

```{r}
set.seed(0)

AR(s = 0.01) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Specify 'along' and 'by' dimensions

```{r, fig.height = 4.5}
set.seed(0)

AR() |>
  generate(n_along = 20, n_by = 3, n_draw = 8) |>
  plot_cor_many()
```

### Specify 'along' and 'by' dimensions, set `con = "by"`

```{r, fig.height = 4.5}
set.seed(0)

AR(con = "by") |>
  generate(n_along = 20, n_by = 3, n_draw = 8) |>
  plot_cor_many()
```


## First-Order Autoregressive `AR1()`

### Model

\begin{equation}
  \beta_j \sim \text{N}\left(\phi \beta_{j-1},  \omega^2\right)
\end{equation}
TMB derives a value of $\omega$ that gives each $\beta_j$ variance $\tau^2$. The prior for $\tau$ is
\begin{equation}
  \tau  \sim \text{N}^+(0, \mathtt{s}^2).
\end{equation}
The prior for $\phi$ is
\begin{equation}
  \frac{\phi - \mathtt{min}}{\mathtt{max} - \mathtt{min}} \sim \text{Beta}(\mathtt{shape1}, \mathtt{shape2}).
\end{equation}

### All defaults

`s = 1`, `min = 0.8`, `max = 0.98`

```{r}
set.seed(0)

AR1() |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Reduce `s`

```{r}
set.seed(0)

AR1(s = 0.01) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Reduce `s` and modify `min`, `max`

```{r}
set.seed(0)

AR1(s = 0.01, min = -1, max = 1) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```



## Random Walk with Seasonal Effects `RW_Seas()`

### Model

\begin{align}
  \beta_j & = \alpha_j + \lambda_j \\
  \alpha_1 & \sim \text{N}(0, \mathtt{sd}^2) \\
  \alpha_j & \sim \text{N}(\alpha_{j-1}, \tau^2), \quad j = 2, \cdots, J \\
  \tau & \sim \text{N}^+\left(0, \mathtt{s}^2\right) \\
  \lambda_j & \sim \text{N}(0, \mathtt{sd\_seas}^2), \quad j = 1, \cdots, \mathtt{n\_seas} - 1 \\
  \lambda_j & = -\sum_{s=1}^{j-1} \lambda_{j-s}, \quad j = \mathtt{n\_seas},\; 2 \mathtt{n\_seas},  \cdots \\
  \lambda_j & \sim \text{N}(\lambda_{j-\mathtt{n\_seas}}, \omega^2), \quad \text{otherwise} \\
  \omega & \sim \text{N}^+\left(0, \mathtt{s\_seas}^2\right)
\end{align}

### All defaults

`s = 1`, `sd = 1`, `s_seas = 0`, `sd_seas = 1` 

```{r}
set.seed(0)

RW_Seas(n_seas = 4) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Reduce `s`, `sd`

`s = 0.01`, `sd = 0`, `s_seas = 0`, `sd_seas = 1` 

```{r}
set.seed(0)

RW_Seas(n_seas = 4, s = 0.01, sd = 0) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Increase `s_seas`

`s = 0.01`, `sd = 0`, `s_seas = 1`, `sd_seas = 1` 

```{r}
set.seed(0)

RW_Seas(n_seas = 4, s = 0.01, sd = 0, s_seas = 1) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Reduce `sd_seas`

`s = 0.01`, `sd = 0`, `s_seas = 1`, `s_seas = 0`, `sd_seas = 0.01` 

```{r}
set.seed(0)

RW_Seas(n_seas = 4, s = 0.01, sd = 0, sd_seas = 0.01) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```



## Second-Order Random Walk with Seasonal Effects `RW2_Seas()`

### Model

\begin{align}
  \beta_j & = \alpha_j + \lambda_j \\
  \alpha_1 & \sim \text{N}(0, \mathtt{sd}^2) \\
  \alpha_2 & \sim \text{N}(\alpha_1, \mathtt{sd\_slope}^2) \\
  \alpha_j & \sim \text{N}(2\alpha_{j-1} - \alpha_{j-2}, \tau^2), \quad j = 3, \cdots, J \\
  \tau & \sim \text{N}^+\left(0, \mathtt{s}^2\right) \\
  \lambda_j & \sim \text{N}(0, \mathtt{sd\_seas}^2), \quad j = 1, \cdots, \mathtt{n\_seas} - 1 \\
  \lambda_j & = -\sum_{s=1}^{j-1} \lambda_{j-s}, \quad j = \mathtt{n\_seas},\; 2 \mathtt{n\_seas},  \cdots \\
  \lambda_j & \sim \text{N}(\lambda_{j-\mathtt{n\_seas}}, \omega^2), \quad \text{otherwise} \\
  \omega & \sim \text{N}^+\left(0, \mathtt{s\_seas}^2\right)
\end{align}

### All defaults

`s = 1`, `sd = 1`, `sd_slope = 1`, `s_seas = 0`, `sd_seas = 1` 

```{r}
set.seed(0)

RW2_Seas(n_seas = 4) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```


### Reduce `s`, `sd`, `sd_slope`, `sd_seas`

`s = 0.01`, `sd = 0`, `sd_slope = 0.01`, `s_seas = 0`, `sd_seas = 0.01` 

```{r}
set.seed(0)

RW2_Seas(n_seas = 4, s = 0.01, sd = 0, sd_slope = 0.01, sd_seas = 0.01) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```


## Linear `Lin()`

### Model

\begin{align}
  \beta_j & = \alpha_j + \epsilon_j \\
  \alpha_j & = \left(j - \frac{J + 1}{2}\right) \eta \\
  \eta & \sim \text{N}\left(\mathtt{mean\_slope}, \mathtt{sd\_slope}^2 \right) \\
  \epsilon & \sim \text{N}(0, \tau^2) \\
  \tau & \sim \text{N}^+\left(0, \mathtt{s}^2\right)
\end{align}

### All defaults

`s = 1`, `mean_slope = 0`, `sd_slope = 1`

```{r}
set.seed(0)

Lin() |>
  generate(n_draw = 8) |>
  plot_cor_one()
```


### Reduce `s`

`s = 0`, `mean_slope = 0`, `sd_slope = 1`

```{r}
set.seed(0)

Lin(s = 0) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Modify `mean_slope`, `sd_slope`

`s = 1`, `mean_slope = 0.2`, `sd_slope = 0.1`

```{r}
set.seed(0)

Lin(mean_slope = 0.2, sd_slope = 0.1) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```


### Specify 'along' and 'by' dimensions

```{r, fig.height = 4.5}
set.seed(0)

Lin() |>
  generate(n_along = 20, n_by = 3, n_draw = 8) |>
  plot_cor_many()
```

### Specify 'along' and 'by' dimensions, set `con = "by"`

```{r, fig.height = 4.5}
set.seed(0)

Lin(con = "by") |>
  generate(n_along = 20, n_by = 3, n_draw = 8) |>
  plot_cor_many()
```


## Linear with AR Errors `Lin_AR()`

### Model

\begin{align}
  \beta_j & = \alpha_j + \epsilon_j \\
  \alpha_j & = \left(j - \frac{J + 1}{2}\right) \eta \\
  \eta & \sim \text{N}\left(\mathtt{mean\_slope}, \mathtt{sd\_slope}^2 \right) \\
  \epsilon_j & \sim \text{N}\left(\phi_1 \epsilon_{j-1} + \cdots + \phi_{\mathtt{n\_coef}} \epsilon_{j-\mathtt{n\_coef}},  \omega^2\right)
  \tau & \sim \text{N}^+\left(0, \mathtt{s}^2\right) \\
  \frac{\phi_k + 1}{2} & \sim \text{Beta}(\mathtt{shape1}, \mathtt{shape2})
\end{align}

### All defaults

`n_coef = 2`, `s = 1`, `mean_slope = 0`, `sd_slope = 1`

```{r}
set.seed(0)

Lin_AR() |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Reduce `s`, `sd_slope`

`s = 0.1`, `mean_slope = 0`, `sd_slope = 0.01`

```{r}
set.seed(0)

Lin_AR(s = 0.1, sd_slope = 0.01) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```


## Linear with AR1 Errors `Lin_AR1()`

### Model

\begin{align}
  \beta_j & = \alpha_j + \epsilon_j \\
  \alpha_j & = \left(j - \frac{J + 1}{2}\right) \eta \\
  \eta & \sim \text{N}\left(\mathtt{mean\_slope}, \mathtt{sd\_slope}^2 \right) \\
  \epsilon_j & \sim \text{N}\left(\phi \epsilon_{j-1},  \omega^2\right)
  \tau & \sim \text{N}^+\left(0, \mathtt{s}^2\right) \\
  \frac{\phi - \mathtt{min}}{\mathtt{max} - \mathtt{min}} \sim \text{Beta}(\mathtt{shape1}, \mathtt{shape2})
\end{align}

### All defaults

`s = 1`, `min = 0.8`, `max = 0.98`, `mean_slope = 0`, `sd_slope = 1`

```{r}
set.seed(0)

Lin_AR1() |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Modify `min`, `max`

`s = 1`, `min = -1`, `max = 1`, `mean_slope = 0`, `sd_slope = 1`

```{r}
set.seed(0)

Lin_AR1(min = -1, max = 1) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```


### Reduce `s`, `sd_slope`

`s = 0.1`, `min = 0.8`, `max = 0.98`, `mean_slope = 0`, `sd_slope = 0.02`

```{r}
set.seed(0)

Lin_AR1(s = 0.1, sd_slope = 0.02) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```


## Penalised Spline `Sp()`

### Model

\begin{align}
  \pmb{\beta} & = \bm{X} \pmb{\alpha} \\
  \alpha_1 & \sim \text{N}(0, \mathtt{sd}^2) \\
  \alpha_2 & \sim \text{N}(\alpha_1, \mathtt{sd\_slope}^2) \\
  \alpha_j & \sim \text{N}(2\alpha_{j-1} - \alpha_{j-2}, \tau^2), \quad j = 3, \cdots, J \\
  \tau & \sim \text{N}^+\left(0, \mathtt{s}^2\right) \\
\end{align}

### All defaults

`n_comp = NULL`, `s = 1`, `sd = 1`, `sd_slope = 1`

```{r}
set.seed(0)

Sp() |>
  generate(n_draw = 8) |>
  plot_cor_one()
```


### Specify `n_comp`

`n_comp = 15`, `s = 1`, `sd = 1`, `sd_slope = 1`

```{r}
set.seed(0)

Sp(n_comp = 5) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

### Reduce `s`, `sd`, `sd_slope`

`n_comp = NULL`, `s = 0.01`, `sd = 0.01`, `sd_slope = 0.01`

```{r}
set.seed(0)

Sp(s = 0.01, sd = 0.01, sd_slope = 0.01) |>
  generate(n_draw = 8) |>
  plot_cor_one()
```

# SVD-Based Priors

## Exchangeable `SVD()`

### Model

\begin{equation}
  \pmb{\beta} = \pmb{F} \pmb{\alpha} + \pmb{g}
\end{equation}

### All defaults

`n_comp = NULL`, `indep = TRUE`

```{r, fig.height = 3}
set.seed(0)

SVD(HMD) |>
  generate(n_draw = 8) |>
  plot_svd_one()
```


### Increase `n_comp`

`n_comp = 5`, `indep = TRUE`

```{r, fig.height = 2}
SVD(HMD, n_comp = 5) |>
  generate(n_draw = 8) |>
  plot_svd_one()
```

### Set `indep` to `FALSE`

```{r, fig.height = 2}
SVD(HMD, indep = FALSE) |>
  generate(n_draw = 8) |>
  plot_svd_one()
```

### Multiple units

```{r, fig.height = 4.5}
SVD(HMD, indep = FALSE) |>
  generate(n_draw = 8, n_element = 3) |>
  plot_svd_many()
```


## Dynamic SVD Prior: `SVD_AR()`

```{r, fig.height = 8}
SVD_AR(HMD, indep = FALSE, s = 0.1) |>
  generate(n_draw = 6, n_along = 5) |>
    ggplot(aes(x = age_mid(age), 
                      y = value, 
                      color = sexgender)) +
      facet_grid(vars(draw), vars(along)) +
      geom_line() +
      scale_color_manual(values = c("darkgreen", "darkorange")) +
      xlab("Age") +
      ylab("") +
      theme(text = element_text(size = 8),
            legend.position = "top",
            legend.title = element_blank())
```