bayestalk.Rmd

---
title: "Introduction to Computational Bayesian Methods for Actuaries"
subtitle: ""
author: "Michael Jones"
date: "10 October 2018"
output:
  xaringan::moon_reader:
    css: "presentation.css"
    lib_dir: libs
    nature:
      highlightStyle: github
      highlightLines: true
      countIncrementalSlides: false
---

```{r setup, include=FALSE}
options(htmltools.dir.version = FALSE)
library(ggplot2)
library(dplyr)
library(svglite)
library(magrittr)
library(tidyr)

knitr::opts_chunk$set(
  dev = "svglite",
  fig.ext = ".svg"
)

nruns <- 1000000 # metropolis example run numbers, increase at final compilation.

```
class: center, inverse, middle


.emph_light[Actuaries]<br>
.emph_dark[should be]<br>
.emph_light[scientific]<br>

---

class: center, inverse, middle


.emph_light[Actuaries]<br>
.emph_dark[should be]<br>
.emph_light[Bayesian]<br>

---

class: center, middle

.emph_blue[What?]

---

class: center, middle

.emph_blue[Why?]

---

class: center, middle

.emph_blue[How?]

---

class: middle, inverse

.emph_dark[Bayesian]<br>
.emph_dark[Recap]

---
class: middle, center

<span style="font-size:600%">
$$P(\theta | D) = \frac{P(D|\theta)P(\theta)}{\int d\theta' P(D|\theta')P(\theta')}$$
</span>

---
class: middle, center

<span style="font-size:600%">
$$\color{lightgrey}{P(\theta | D) = \frac{P(D|\theta)\color{black}{P(\theta)}}{\int d\theta' P(D|\theta')P(\theta')}}$$
</span>
---
class: middle, center

<span style="font-size:600%">
$$\color{lightgrey}{P(\theta | D) = \frac{\color{black}{P(D|\theta)}{P(\theta)}}{\int d\theta' P(D|\theta')P(\theta')}}$$
</span>

---

class: middle, center

<span style="font-size:600%">
$$\color{lightgrey}{\color{black}{P(\theta | D)} = \frac{P(D|\theta){P(\theta)}}{\int d\theta' P(D|\theta')P(\theta')}}$$
</span>

---

class: middle, center

<span style="font-size:600%">
$$\color{lightgrey}{P(\theta | D) = \frac{P(D|\theta){P(\theta)}}{\color{black}{\int d\theta' P(D|\theta')P(\theta')}}}$$
</span>

---



class: middle, inverse

.emph_dark[Why be]<br>
.emph_dark[Bayesian]


---

class: middle

.emph_blue[It works]

???

It produces results that are testable and given the right models generally are as good (sometimes better) than frequentist analysis.

---

class: middle

.emph_blue[It's coherent]

???

Most Bayesian analysis stems from a few foundational ideas compared to frequentist statistics which has a much more varied set of fundamentals.

Posterior information is rich and useful

---
class: middle

.emph_blue[It's natural]

???

Humans naturally think in Bayesian terms which are then 'unlearned' when studying frequentist analysis. 


---

class: middle

.emph_blue[It's self]<br>
.emph_blue[documenting]

???

Through the priors, assumptions are formally absorbed into the analysis and cannot be hidden as in classical statistics.

---

class: middle

.emph_blue[You get the]<br>
.emph_blue[full Posterior]


---

class: middle

.emph_blue[It's modular]


???

Multilevel models, choice of priors, e.g. simple to make robust to outliers


---

class: middle, inverse

.emph_dark[Why ] .emph_light[not ] .emph_dark[be]<br>
.emph_dark[Bayesian]

---
class: middle

.emph_blue[What]<br>
.emph_blue[prior?]

---

class: middle

.emph_blue[It's]<br>
.emph_blue[awkward]

???

A whole posterior distribution can be difficult to communicate/process

---

class: middle

.emph_blue[It's]<br>
.emph_blue[hard]

???

ok if conjugate pairs, but otherwise, it's difficult to do without MCMC or other computationally intensive programs.


---

class: middle

.emph_blue[It's]<br>
.emph_blue[weird]

???

It's not what people are used to, and the themes are unfamiliar, though it's gaining traction in social sciences and pharmacology (and election prediction)

---

class: middle, inverse

.emph_light[Frequentism]<br>
.emph_dark[vs]<br>
.emph_light[Bayesianism]
---

class: middle

.emph_blue[Probability]

---
class: middle

.emph_blue[Data]


---
class: middle

.emph_blue[Parameters]


---

class: middle

.emph_blue[Ranges]

---


```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height=2}
# From https://freakonometrics.hypotheses.org/18117
set.seed(200)
xbar <- 150
n <- 1000
ns <- 100
M <- matrix(rbinom(n*ns, size=1, prob=xbar/n),
            nrow=n)

fIC <- function(x) {
  mean(x) + c(-1,1)*1.96*sqrt(mean(x)*(1-mean(x)))/sqrt(n)
}

IC <- t(apply(M,2,fIC))
MN <- apply(M,2,mean)

k <- (xbar/n<IC[,1])|(xbar/n>IC[,2])
#plot(MN, 1:ns, xlim=range(IC), axes=FALSE,
#     xlab='', ylab='', pch=19, cex=.7,
#     col=c('blue','red')[1+k])
#axis(1)
#segments(IC[,1],1:ns,IC[,2],1:ns,
#         col=c('blue','red')[1+k])
#abline(v=xbar/n)

to_plot <- data.frame(MN, left=IC[,1], right=IC[,2],k, y=1:ns)

s1 <- svgstring(standalone = FALSE)
ggplot(to_plot) +
  geom_point(aes(MN, y, col=k)) +
  geom_errorbarh(aes(xmin=left, xmax=right, y=y, colour=k), height=0) +
  geom_vline(xintercept=xbar/n) +
  scale_y_continuous('') +
  ggtitle('Frequentist Confidence Intervals', subtitle='100 draws from a Binomial with mean 15') +
  scale_x_continuous('') +
  theme_minimal() +
  theme(legend.position = 'none',
        axis.text.y = element_blank())

htmltools::HTML(s1())
invisible(dev.off())
```

---

```{r echo=FALSE, message=FALSE, warning=FALSE}


y <- rbeta(10000, 3,7)
bayes_range <- data.frame(y=y)

hdi <- coda::HPDinterval(coda::as.mcmc(y), prob=0.95)
label_x_pos <- hdi[1] + 0.75*(hdi[2]-hdi[1])
label_y_pos <- 760
label_text <- paste('95% HDI:', round(hdi[1],2), 'to', round(hdi[2],2), sep=' ')

label_df <- data.frame(x=label_x_pos, y=label_y_pos, text=label_text)


s <- svgstring(standalone=FALSE)
ggplot(bayes_range) + 
  geom_histogram(aes(y), fill='#30499b', colour='white') +
  geom_vline(xintercept = hdi[1], linetype='dashed') +
  geom_vline(xintercept = hdi[2], linetype='dashed') +
  geom_label(aes(x=x, y=y, label=text), data=label_df) +
  ggtitle('Bayesian Highest Density Intervals') +
  theme_bw()
htmltools::HTML(s())

```



---

class: middle

.emph_blue[p-values &]<br>
.emph_blue[significance]<br>
.emph_blue[testing]



---

class: middle, inverse

.emph_dark[How to be]<br>
.emph_dark[Bayesian]

---

class: middle

.emph_blue[The Bayesian]<br>
.emph_blue[Workflow]

---

class: middle

- Identify your data
- Define a descriptive model
- Specify a prior
- Compute the Posterior
- Interpret the Posterior
- Check the model is reasonable

---

class: middle

.emph_dark[Analytically:]<br>
.emph_blue[Estimating the]<br>
.emph_blue[Bias of a coin]

---

##Parameter

<span style="font-size:400%">
$$\theta = P(y_i = \text{heads})$$
</span>

---

## Likelihood


### Bernoulli:

<span style="font-size:400%">
$$ p(y|\theta) = \theta^y(1 - \theta)^{(1-y)}$$
</span>

---

## Prior
### Beta
<span style="font-size:400%">
$$\begin{align}
  P(\theta|a,b) &= \text{beta}(a,b)\\\\
  &=\frac{\theta^{(a-1)}(1-\theta)^{(b-1)}}{B(a,b)}
  \end{align}$$
</span>

---


```{r echo=FALSE, message=FALSE, warning=FALSE, fig.width=12, fig.height=9, dpi=600}
a_vals <- c(0.1,1,2,3,4,20)
b_vals <- c(0.1,1,2,3,4,20)
theta_vals <- seq(0,1,0.01)

theta <- rep(theta_vals, length(a_vals) * length(b_vals))
a <- rep(a_vals, each = length(b_vals) * length(theta_vals))
b <- rep(rep(b_vals, each = length(theta_vals)), length(a_vals))

plot_data <- data.frame(a,b,theta)
plot_data <- mutate(plot_data, y=dbeta(theta, a, b))


s <- svgstring(standalone=FALSE)
ggplot(plot_data, aes(theta, y)) +
  geom_line() +
  facet_wrap(a~b, scales='free_y',labeller=label_wrap_gen(multi_line=FALSE)) +
  scale_x_continuous('θ', breaks=c(0,0.25,0.5,0.75,1)) + 
  theme_bw() +
  theme(axis.text=element_blank(),
        axis.ticks=element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank())

htmltools::HTML(s())


```


---
## Posterior
### Another Beta

<span style="font-size:300%">
$$\begin{align}
  P(\theta|n_{heads},n_{tails},a,b) &= \text{beta}(a+n_{heads},b+n_{tails})\\\\
  &=\frac{\theta^{(a+n_{heads}-1)}(1-\theta)^{(b+n_{tails}-1)}}{B(a+n_{heads},b+n_{tails})}
  \end{align}$$
</span>

---

# Live demo time:

https://mjones.shinyapps.io/coin/

.pull-left[
## If it works

click [here](#47)
]

.pull-right[
## If it doesn't work

click [here](#43)
]

---

# No idea about the coin
Flat $\text{beta}(1,1)$ distribution:

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height=1, fig.width=2}
theta <- seq(0,1,0.1)
df <- data.frame(theta, y=dbeta(theta, 1, 1))

s <- svgstring(standalone=FALSE, height=5)
ggplot(df, aes(theta, y)) + 
  geom_line() +
  theme_bw() +
  scale_x_continuous('θ') +
  scale_y_continuous('Prior')
htmltools::HTML(s())
```

---

# Collect some data

7 Heads, 5 tails:

```{r echo=FALSE}
theta <- seq(0,1,0.01)
likelihood <- theta^7 * (1 - theta) ^5

s <- svgstring(standalone = FALSE, height=5)
df <- data.frame(theta, likelihood)
ggplot(df, aes(theta, likelihood)) +
  geom_line() +
  theme_bw()
htmltools::HTML(s())

```

---
# Posterior

Another Beta: $Beta(6,8)$

```{r echo=FALSE, message=FALSE, warning=FALSE}
theta <- seq(0,1,0.01)
posterior <- theta^8 * (1 - theta) ^6

s <- svgstring(standalone = FALSE, height=5)
df <- data.frame(theta, posterior)
ggplot(df, aes(theta, posterior)) +
  geom_line() +
  theme_bw()
htmltools::HTML(s())
```
---

# Strong Prior

```{r echo=FALSE}
a <- 2
b <- 20
n_heads <- 7
n_tails <- 5

data <- data.frame(theta=seq(0,1,0.001))
data %<>% mutate(prior=dbeta(theta, shape1=a, shape2=b),
                 likelihood=theta^n_heads * (1 - theta) ^n_tails,
                 posterior=dbeta(theta, shape1=a+n_heads, shape2=b+n_tails))

data %<>% gather(type,y, prior, likelihood, posterior) %>%
  mutate(type = ordered(type, levels=c('prior','likelihood','posterior')))

s <- svgstring(standalone=FALSE, height=5)
ggplot(data, aes(theta, y)) +
  geom_line() +
  facet_grid(type~., scales='free_y') +
  theme_bw()
htmltools::HTML(s())
```

---

class: middle

.emph_blue[Problems]

???

- Conjugate priors
- calculating the normalisation integral

---

class: middle

.emph_blue[What prior]<br>
.emph_blue[to use?]

---

class: middle

.emph_blue[Normalising]

---

class: middle

.emph_dark[Computationally:]<br>
.emph_blue[First Steps]
---

```{r echo=FALSE}
x <- seq(0,1,0.01)
grid_beta_a <- 3
grid_beta_b <- 6
y <- dbeta(x, grid_beta_a, grid_beta_b)
grid_plot_df <- data.frame(x,y)

grid_plot <- ggplot() + geom_line(data=grid_plot_df, aes(x,y)) + theme_bw()

s <- svgstring(standalone = FALSE)
grid_plot
htmltools::HTML(s())

```
---

```{r echo=FALSE}
x_grid <- seq(0,1,0.1)
y_grid_base <- rep(0,length(x_grid))

grid_base_df <- data.frame(x_grid, y_grid_base)

grid_plot <- grid_plot + geom_point(data=grid_base_df, aes(x_grid, y_grid_base))

s <- svgstring(standalone = FALSE)
grid_plot
htmltools::HTML(s())
```
---
```{r echo=FALSE}
y_grid_tops <- dbeta(x_grid, grid_beta_a, grid_beta_b) / 1.4

grid_base_df$y_grid_tops <- y_grid_tops

grid_plot <- ggplot() + geom_line(data=grid_plot_df, aes(x,y), colour='grey') +
  geom_point(data=grid_base_df, aes(x_grid, y_grid_base)) +
  geom_point(data=grid_base_df, aes(x_grid, y_grid_tops)) +
  geom_segment(data = grid_base_df, aes(x=x_grid, y=y_grid_base, xend=x_grid, yend=y_grid_tops)) +
  theme_bw()

s <- svgstring(standalone = FALSE)
grid_plot
htmltools::HTML(s())
```
---
```{r echo=FALSE}
y_grid_tops <- dbeta(x_grid, grid_beta_a, grid_beta_b)

grid_base_df$y_grid_tops <- y_grid_tops

grid_plot <- ggplot() + geom_line(data=grid_plot_df, aes(x,y), colour='grey') +
  geom_point(data=grid_base_df, aes(x_grid, y_grid_base)) +
  geom_point(data=grid_base_df, aes(x_grid, y_grid_tops)) +
  geom_segment(data = grid_base_df, aes(x=x_grid, y=y_grid_base, xend=x_grid, yend=y_grid_tops)) +
  theme_bw()

s <- svgstring(standalone = FALSE)
grid_plot
htmltools::HTML(s())

```


---

```{r echo=FALSE}
x <- seq(0,1,0.001)
grid_beta_a <- 300
grid_beta_b <- 600
y <- dbeta(x, grid_beta_a, grid_beta_b)
grid_plot_df <- data.frame(x,y)

grid_plot <- ggplot() + geom_line(data=grid_plot_df, aes(x,y)) + 
  geom_point(data=grid_base_df, aes(x_grid, y_grid_base))+
  theme_bw()

s <- svgstring(standalone = FALSE)
grid_plot
htmltools::HTML(s())

```

???

- If the posterior is narrower, you waste a lot of time/computing power calculating the posterior in places which do not matter
- Not only that, but as your dimensions increase (i.e. number of parameters), the number of points increases exponentially.

---


class: middle

.emph_dark[Computationally:]<br>
.emph_blue[MCMC]

---

class:middle

.emph_blue[Markov Chain]<br>
.emph_blue[Monte Carlo]
---

class: middle

.emph_dark[Markov Chain]<br>
.emph_blue[Monte Carlo]

---

class: middle

.emph_blue[Markov Chain]<br>
.emph_dark[Monte Carlo]

---

class: middle

.emph_blue[Metropolis]


---
# Metropolis Algorithm
```{r echo=FALSE}
knitr::include_graphics('metropolis_pics/1_current_place.svg')
```

---
# Metropolis Algorithm
```{r echo=FALSE}
knitr::include_graphics('metropolis_pics/2_choices.svg')
```

---
# Metropolis Algorithm
```{r echo=FALSE}
knitr::include_graphics('metropolis_pics/3_chosen.svg')
```

---
# Metropolis Algorithm
```{r echo=FALSE}
knitr::include_graphics('metropolis_pics/4_definitely_move.svg')
```

---
# Metropolis Algorithm
```{r echo=FALSE}
knitr::include_graphics('metropolis_pics/5_maybe_move.svg')
```
---
# Metropolis Algorithm Example

--

- 50 Islands

--

- Populations in ratio 1:2:3:...:50

--

- Want to visit each in accordance with its population



---

```{r metropolis_sim, include=FALSE}
n_islands <- 50
islands <- c(0,1:n_islands,0)
island_density <- islands / sum(islands)

starting_pos <- 2
print(starting_pos)

visited = vector(length=nruns)

decide_candidate <- function(current_pos) {
  if (runif(1) < 0.5) {
    return(current_pos - 1)
  } else {
    return(current_pos + 1)
  }
}


decide_to_move <- function(current_pos, candidate) {
  #print(paste('current position: ', current_pos, ' population: ', islands[current_pos], sep=''))
  #print(paste('candidate: ', candidate, ' population: ', islands[candidate], sep=''))
  if (islands[candidate] > islands[current_pos]) {
    return(TRUE)
  } else {
    if (runif(1) < islands[candidate] / islands[current_pos]) {
      return(TRUE)
    } else {
      return(FALSE)
    }
  }
}

jump <- function(current_pos) {
  candidate <- decide_candidate(current_pos)
  
  if (decide_to_move(current_pos, candidate)) {
    return(candidate)
  } else {
    return(current_pos)
  }
}

# run the sim
# initialise
visited[1] <- starting_pos


for (i in 2:nruns) {
  visited[i] = jump(visited[i-1])
}

# re label the islands 1:7
visited <- visited - 1

freq_table <- function(visited) {
  #df <- as.data.frame(table(visited))
  df <- data.frame(visited = 2:(length(islands) - 1))
  
  Freq <- vector(length=length(df$visited))
  for (i in 1:length(Freq)) {
    Freq[i] <- sum(visited == df$visited[i])
  }
  df$Freq <- Freq
  
  df <- mutate(df, density = Freq/sum(Freq))
  
  return(df)
}


```

```{r echo=FALSE}
df <- data.frame(index=1:length(visited), visited=visited)


#ggplot(df, aes(index, visited)) + geom_line() +  scale_x_log10()
df %<>% mutate(slice = ifelse(index < 500, 'beginning', ifelse(index > max(index) - 500, 'end','other'))) %>% filter(slice != 'other')

s <- svgstring(standalone = FALSE)
ggplot(df, aes(index, visited)) +
  geom_line() +
  facet_wrap(~slice, scales='free_x') +
  scale_y_continuous('Island Number') +
  scale_x_continuous('Step Number') +
  theme_bw() 
htmltools::HTML(s())
```

---

```{r echo=FALSE}
frequency_matrix <- function(visited, stages, candidates) {
  results <- matrix(nrow=length(candidates), ncol=length(stages))
  for (i in seq_along(stages)) {
    results[,i] <-freq_dist(visited[1:stages[i]], candidates)
  }
  
  colnames(results) <- stages
  results <- as.data.frame(results)
  results$island <- candidates
  results <- gather(results, stage, value, -island)
  results$stage <- factor(results$stage, levels=ordered(stages))
  return(results)
}


freq_dist <- function(visited, candidates) {
  frequencies <- vector(length=length(candidates))
  for (i in seq_along(candidates)) {
    frequencies[i] <- sum(visited == candidates[i])
  }
  return(frequencies/sum(frequencies))
}

candidates <- 1:n_islands
stages <- c(1, 10, 20, 100, 200, 300, 1000, 10000, nruns)
islands_frequencies <- frequency_matrix(visited, stages, candidates)

s <- svgstring(standalone = FALSE)
ggplot(islands_frequencies) +
  geom_col(aes(island, value)) +
  facet_wrap(~stage, scales='free_y') +
  scale_y_continuous('Visiting proportion') +
  scale_x_continuous('Step Number') +
  theme_bw()
htmltools::HTML(s())
```


---

class: middle

.emph_blue[Don't make]<br>
.emph_blue[your own]


---

class: middle

.emph_dark[Computationally:]<br>
.emph_blue[In practice]

---
# Two Coins
```{r echo=FALSE}
knitr::include_graphics('two_coins.svg')
```

---
class: center, middle 
```{r message=FALSE, warning=FALSE, include=FALSE}
s <- rep(c(1,2), each=10)
y <- c(1,1,1,1,1,1,1,0,0,0,1,1,1,1,0,0,0,0,0,0)

two_coin_data <- data.frame(s,y)

source('../DBDA/Jags-Ydich-XnomSsubj-MbernBeta.R')
mcmcCoda <- genMCMC(data=two_coin_data, numSavedSteps=10000)
```


```{r echo=FALSE}
head(mcmcCoda[[1]], n=30)
```

---
```{r echo=FALSE, fig.height=6}
s <- svgstring(standalone = FALSE)
plotMCMC(mcmcCoda, data=two_coin_data, compVal = NULL, compValDiff = 0.0)
htmltools::HTML(s())
```


---

class: middle, inverse

.emph_dark[Bayesianism]<br>
.emph_dark[applied to]<br>
.emph_dark[Insurance]

---

class: middle, center

```{r echo=FALSE}
dat <- read.csv('ClarkTriangle.csv')
options(knitr.kable.NA = '')
triangle_table <- select(dat, -X) %>% spread(dev, cum) %>% mutate(AY = factor(AY))
knitr::kable(triangle_table, format = 'html', format.args=list(big.mark=','), digits=0)
```
---
```{r echo=FALSE}
s <- svgstring(standalone = FALSE)
ggplot(dat, aes(dev, cum, colour=factor(AY))) + geom_point() + geom_line() + theme_bw() + scale_x_continuous('Development Period') + scale_y_continuous('Amount') + theme(legend.title=element_blank())
htmltools::HTML(s())
```

---
# The model
See [here](https://magesblog.com/post/2015-11-10-hierarchical-loss-reserving-with-stan/)
$$\begin{align}
CL_{AY,dev} &\sim \mathcal{N}(\mu_{AY, dev}, \sigma^2_{dev})\\
\mu_{AY,dev} &= ULT_{AY} \cdot G(dev|\omega, \theta)\\
\sigma_{dev}&=\sigma \sqrt{\mu_{dev}}\\
ULT_{AY} &\sim \mathcal{N}(\mu_{ult}, \sigma^2_{ult})\\
G(dev|\omega,\theta) &=1-\exp\left(-\left(\frac{dev}{\theta}\right)^{\omega}\right)
\end{align}$$
---
```{r include=FALSE}
dat$origin <- dat$AY-min(dat$AY)+1
dat <- dat[order(dat$dev),]
# Add future dev years
nyears <- 12
newdat <- data.frame(
  origin=rep(1:10, each=nyears),
  AY=rep(sort(unique(dat$AY)), each=nyears),  
  dev=rep(seq(from=6, to=nyears*12-6, by=12), 10)
)
newdat <- merge(dat, newdat, all=TRUE)
newdat <- newdat[order(newdat$dev),]
start.vals <- c(ult = 5000, omega = 1.4, theta = 45)

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores()) 
mlstanDso <- stan_model(file = "Hierachical_Model_for_Loss_Reserving.stan",
                        model_name = "MultiLevelGrowthCurve") 
newdat <- newdat[is.na(newdat$cum),]

mlfit <- sampling(mlstanDso, iter=7000, warmup=2000, 
                  thin=2, chains=4, 
                  data = list(N = nrow(dat),
                              cum=dat$cum,
                              origin=dat$origin,
                              n_origin=length(unique(dat$origin)),
                              dev=dat$dev,
                              N_mis=nrow(newdat),
                              origin_mis=newdat$origin,
                              dev_mis=newdat$dev))
```

# Stan model output

```{r echo=FALSE, fig.height=10}
print(mlfit, c("mu_ult", "omega", "theta", "sigma_ult", "sigma"),
probs=c(0.5, 0.75, 0.975))
```

---


```{r echo=FALSE}
ult <- as.data.frame(rstan::extract(mlfit, paste0("ult[", 1:10, "]")))

ult <- gather(ult, year, value)
ult$year <- factor(ult$year, levels=c("ult.1.","ult.2.","ult.3.", "ult.4.", "ult.5.", "ult.6.", "ult.7.", "ult.8.", "ult.9.", "ult.10."))

s <- svgstring(standalone = FALSE)
ggplot(ult) + geom_density(aes(value), fill='lightblue') + facet_wrap(~year, scale='free') + theme_bw()
htmltools::HTML(s())
```

---

```{r echo=FALSE, message=FALSE, warning=FALSE}
Y_mean <- rstan::extract(mlfit, "Y_mean")
Y_mean_cred <- apply(Y_mean$Y_mean, 2, quantile, c(0.025, 0.975)) 
Y_mean_mean <- apply(Y_mean$Y_mean, 2, mean)
Y_pred <- rstan::extract(mlfit, "Y_pred")
Y_pred_cred <- apply(Y_pred$Y_pred, 2, quantile, c(0.025, 0.975)) 
Y_pred_mean <- apply(Y_pred$Y_pred, 2, mean)
dat2 <- rbind(dat, newdat)
dat2$Y_pred_mean <- Y_pred_mean
dat2$Y_pred_cred5 <- Y_pred_cred[1,]
dat2$Y_pred_cred95 <- Y_pred_cred[2,]

s <- svgstring(standalone = FALSE)
ggplot(dat2) +
  geom_ribbon(aes(x=dev, ymin=Y_pred_cred5, ymax=Y_pred_cred95), fill='lightblue', alpha=0.5) +
  geom_point(aes(dev, cum)) +
  geom_line(aes(dev, Y_pred_mean)) + 
  facet_wrap(~factor(AY)) +
  theme_bw()
htmltools::HTML(s())
```
---

# Shinystan
Available [here](https://mjones.shinyapps.io/staninsurance/)


---

class: middle inverse

.emph_dark[Further]<br>
.emph_dark[Topics]


---

class: middle

.emph_blue[Multilevel]<br>
.emph_blue[Models]

---

class: middle

.emph_blue[Bayesian]<br>
.emph_blue[Model]<br>
.emph_blue[Averaging]

---

class: middle

.emph_blue[Causal]<br>
.emph_blue[Networks]

---

# References and Further Reading

- *Doing Bayesian Data Analysis* by John K. Kruschke
- *Computational Actuarial Science* edited by Arthur Charpentier
- *Bayesian Data Analysis (3rd edition)* by Andrew Gelman *et al*
- Markus Gesmann's Blog (https://magesblog.com/)
- Arthur Charpentier's blog (https://freakonometrics.hypotheses.org/)

# Interactive

- http://mjones.shinyapps.io/coin/
- http://mjones.shinyapps.io/staninsurance/