Practical 8

class: center, middle, title-slide

.title[
# Practical 8
]
.author[
### Olivier Gimenez
]
.date[
### last updated: 2024-10-17
]

---

## Model selection with wAIC

* Consider the plant example.

* Compare the three models (no, partial and complete pooling) with wAIC.

---

# Solution

---
## wAIC of model with no pooling

``` r
samples <- jags.samples(model = allom.1$model,
                        variable.names = c("WAIC", "deviance"), 
                        type = "mean", 
                        n.iter = 2000, 
                        n.burnin = 1000, 
                        n.thin = 1)
samples$p_waic <- samples$WAIC
samples$waic <- samples$deviance + samples$p_waic
tmp <- sapply(samples, sum)
waic_completepooling <- round(c(waic = tmp[["waic"]], 
                                p_waic = tmp[["p_waic"]]),1)
```

---
## wAIC of model with partial pooling

``` r
samples <- jags.samples(model = allom.2$model,
                        variable.names = c("WAIC","deviance"), 
                        type = "mean", 
                        n.iter = 2000, 
                        n.burnin = 1000, 
                        n.thin = 1)
samples$p_waic <- samples$WAIC
samples$waic <- samples$deviance + samples$p_waic
tmp <- sapply(samples, sum)
waic_partialpooling <- round(c(waic = tmp[["waic"]], 
                               p_waic = tmp[["p_waic"]]),1)
```

---
## wAIC of model with complete pooling

``` r
samples <- jags.samples(model = allom.3$model,
                        variable.names = c("WAIC","deviance"), 
                        type = "mean", 
                        n.iter = 2000, 
                        n.burnin = 1000, 
                        n.thin = 1)
samples$p_waic <- samples$WAIC
samples$waic <- samples$deviance + samples$p_waic
tmp <- sapply(samples, sum)
waic_nopooling <- round(c(waic = tmp[["waic"]], 
                          p_waic = tmp[["p_waic"]]),1)
```

---
## Model ranking

``` r
data.frame(model = c('no pooling', 
                     'partial pooling', 
                     'complete pooling'),
           waic = c(waic_nopooling[1],
                    waic_partialpooling[1],
                    waic_completepooling[1]),
           p_waic = c(waic_nopooling[2],
                      waic_partialpooling[2],
                      waic_completepooling[2])) %>%
  arrange(waic)
```

---
## Model ranking

```
##              model   waic p_waic
## 1       no pooling 2405.8   65.8
## 2  partial pooling 2524.4   45.9
## 3 complete pooling 3674.4    2.4
```