bglmm_example4.qmd

---
title: "Bayesian GLMM Part4"
author: "Murray Logan"
date: today
date-format: "DD/MM/YYYY"
format: 
  html:
    ## Format
    theme: [default, ../resources/ws-style.scss]
    css: ../resources/ws_style.css
    html-math-method: mathjax
    ## Table of contents
    toc: true
    toc-float: true
    ## Numbering
    number-sections: true
    number-depth: 3
    ## Layout
    page-layout: full
    fig-caption-location: "bottom"
    fig-align: "center"
    fig-width: 4
    fig-height: 4
    fig-dpi: 72
    tbl-cap-location: top
    ## Code
    code-fold: false
    code-tools: true
    code-summary: "Show the code"
    code-line-numbers: true
    code-block-border-left: "#ccc"
    code-copy: true
    highlight-style: atom-one
    ## Execution
    execute:
      echo: true
      cache: true
    ## Rendering
    embed-resources: true
crossref:
  fig-title: '**Figure**'
  fig-labels: arabic
  tbl-title: '**Table**'
  tbl-labels: arabic
engine: knitr
output_dir: "docs"
documentclass: article
fontsize: 12pt
mainfont: Arial
mathfont: LiberationMono
monofont: DejaVu Sans Mono
classoption: a4paper
bibliography: ../resources/references.bib
---

```{r}
#| label: setup
#| include: false
#| cache: false

knitr::opts_chunk$set(cache.lazy = FALSE,
                      tidy = "styler")
options(tinytex.engine = "xelatex")
```

# Preparations

Load the necessary libraries

```{r}
#| label: libraries
#| output: false
#| eval: true
#| warning: false
#| message: false
#| cache: false
library(tidyverse) #for data wrangling
library(car)       #for regression diagnostics
library(broom)     #for tidy output
#library(ggfortify) #for model diagnostics
library(knitr)     #for kable
library(emmeans)   #for estimating marginal means
library(MASS)      #for glm.nb
library(brms)
library(broom.mixed)
library(tidybayes)
library(bayesplot)
library(standist)   #for visualizing distributions
library(rstanarm)
library(cmdstanr)   #for cmdstan
library(ggeffects)
library(rstan)
library(DHARMa)
library(ggridges)
library(easystats)     #framework for stats, modelling and visualisation
library(patchwork)
library(modelsummary)
source('helperFunctions.R')
```

# Scenario

![Crab_shrimp_coral](../resources/crab_shrimp_coral.jpg){#fig-crabs width="400" height="284"}

To investigate synergistic coral defence by mutualist crustaceans,
@Mckeon-2012-1095 conducted an aquaria experiment in which colonies of a coral
species were placed in a tank along with a predatory sea star and one of four
symbiont combinations:

- no symbiont,
- a crab symbiont
- a shrimp symbiont
- both a crab and shrimp symbiont.

The experiments were conducted in a large octagonal flow-through seawater tank
that was partitioned into eight sections, which thereby permitted two of each of
the four symbiont combinations to be observed concurrently. The tank was left
overnight and in the morning, the presence of feeding scars on each coral colony
was scored as evidence of predation.  The experiments were repeated ten times,
each time with fresh coral colonies, sea stars and symbiont.

The ten experimental times represent blocks (random effects) within which the
symbiont type (fixed effect) are nested.

# Read in the data

```{r readData, results='markdown', eval=TRUE}
mckeon <- read_csv("../data/mckeon.csv", trim_ws = TRUE)
```


# Exploratory data analysis

```{r}
# Data is presence/absence so the model is binomial, it has logit link. We have intercept, bunch of slopes from different treatment types and another random effect called block (think of it as CoralID or tank)
```

```{r}
mckeon <- mckeon |> mutate(BLOCK = factor(BLOCK), 
SYMBIONT = factor(SYMBIONT, levels = c('none','crabs','shrimp','both')))


mckeon |> 
ggplot(aes(y = PREDATION, x = SYMBIONT)) +
  geom_point(position = position_jitter(width = 0.2, height = 0)) + facet_grid(~BLOCK) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
  
  
# The blocks are all acting very different, we might need to add
```

Model formula:
$$
y_i \sim{} \mathcal{N}(n, p_i)\\
ln\left(\frac{p_i}{1-p_1}\right) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i}
$$

where $\boldsymbol{\beta}$ and $\boldsymbol{\gamma}$ are vectors of the fixed and random effects parameters respectively 
and $\bf{X}$ is the model matrix representing the overall intercept and effects of symbionts on the probability of the colony experiencing predation.
$\bf{Z}$ represents a cell means model matrix for the random intercepts associated with individual coral colonies.


# Fit the model

```{r}
#| label: Formula

mckeon.form <- bf(PREDATION | trials(1) ~ SYMBIONT + (1|BLOCK),
                  family = binomial(link = 'logit'))

```



```{r}
#| label: Get the Priors
get_prior(mckeon.form, data = mckeon)



priors <- prior(normal(0,1.8), class = "Intercept") + 
  prior(normal(0,1), class = "b") +
  prior(student_t(3,0,1.8), class = "sd")
```

```{r}
#| label: Run the chains
mckeon.brm2 <- brm(mckeon.form,
                 data = mckeon,
                 prior = priors,
                 sample_prior = 'only',
                 iter = 5000,
                 warmup = 1000,
                 chains = 3, cores = 3,
                 thin = 5,
                 control = list(adapt_delta = 0.99, max_treedepth = 20),
                 refresh = 0,
                 backend = "cmdstanr")
```

```{r}
mckeon.brm2 |> 
  conditional_effects() |> 
  plot(points = TRUE)
```


```{r}
#| label: Introduce data
mckeon.brm2 <- brm(mckeon.form,
                 data = mckeon,
                 prior = priors,
                 sample_prior = 'yes',
                 iter = 5000,
                 warmup = 1000,
                 chains = 3, cores = 3,
                 thin = 5,
                 control = list(adapt_delta = 0.99, max_treedepth = 20),
                 refresh = 0,
                 backend = "cmdstanr")
```

```{r}
mckeon.brm2 |> 
  conditional_effects() |> 
  plot(points = TRUE)

mckeon.brm3 <- update(mckeon.brm2,sample_prior = 'yes',
                      refresh = 0)
```

```{r}
mckeon.brm3$fit |> stan_trace() #pars not specified can be removed by specifying
mckeon.brm3$fit |> stan_ac()
mckeon.brm3$fit |> stan_rhat()
mckeon.brm3$fit |> stan_ess() 
mckeon.brm3 |>  pp_check(type = 'dens_overlay',, ndraws = 250)

mckeon.resids <- make_brms_dharma_res(mckeon.brm3, integerResponse = FALSE)
testUniformity(mckeon.resids)
plotResiduals(mckeon.resids, quantreg = FALSE)
plot(mckeon.resids, form = factor(rep(1, nrow(mckeon))))
testDispersion(mckeon.resids)
plotResiduals(mckeon.resids, quantreg = TRUE)

mckeon.form <- bf(PREDATION | trials(1) ~ SYMBIONT + (SYMBIONT|BLOCK),
                  family = binomial(link = 'logit'))

priors <- prior(normal(0,1.8), class = "Intercept") + 
  prior(normal(0,3), class = "b") +
  prior(student_t(3,0,3), class = "sd") +
  prior(lkj_corr_cholesky(1), class = 'cor')
```

```{r}
#| label: New model with added prior
mckeon.brm4 <- brm(mckeon.form,
                   data = mckeon,
                   prior = priors,
                   sample_prior = 'yes',
                   iter = 5000,
                   warmup = 2500,
                   chains = 3, cores = 3,
                   thin = 5,
                   refresh = 0,
                   control = list(adapt_delta = 0.99, max_treedepth = 20),
                   backend = 'cmdstanr')

mckeon.brm4 |>  SUYR_prior_and_posterior()

```

# MCMC sampling diagnostics 
```{r}
mckeon.resids <- make_brms_dharma_res(mckeon.brm4, integerResponse = FALSE)
testUniformity(mckeon.resids)
plotResiduals(mckeon.resids, qunatreg = FALSE)
testDispersion(mckeon.resids)

mckeon.brm4 |> loo()
mckeon.brm3 |> loo()
loo_compare(loo(mckeon.brm3), loo(mckeon.brm4))
```



# Model validation 
::: {.panel-tabset}
## brms 
:::: {.panel-tabset}
### pp check
Post predictive checks provide additional diagnostics about the fit of the
model.  Specifically, they provide a comparison between predictions drawn from
the model and the observed data used to train the model.


- dens_overlay: plots the density distribution of the observed data (black line)
overlayed on top of 50 density distributions generated from draws from the model
(light blue).  Ideally, the 50 realisations should be roughly consistent with
the observed data.

The model draws appear to be consistent with the observed data.

- error_scatter_avg: this plots the observed values against the average
  residuals. Similar to a residual plot, we do not want to see any patterns in
  this plot.  Note, this is not really that useful for models that involve a
  binomial response


This is not really interpretable

- intervals:  plots the observed data overlayed on top of posterior predictions
associated with each level of the predictor.  Ideally, the observed data should
all fall within the predictive intervals.



The `shinystan` package allows the full suite of MCMC diagnostics and posterior
predictive checks to be accessed via a web interface.


### DHARMa residuals

DHARMa residuals provide very useful diagnostics.  Unfortunately, we cannot
directly use the `simulateResiduals()` function to generate the simulated
residuals.  However, if we are willing to calculate some of the components
yourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

- simulated (predicted) responses associated with each observation.
- observed values
- fitted (predicted) responses (averaged) associated with each observation


**Conclusions:**

- the simulated residuals do not suggest any issues with the residuals
- there is no evidence of a lack of fit.

::::
:::
<!-- END_PRIVATE-->
# Partial effects plots 

```{r}
mckeon.brm4 |> 
  conditional_effects() |> 
  plot(points = TRUE)

mckeon.brm4 |> 
  as_draws_df() |> 
  colnames()
```


```{r}
mckeon.brm4 |> 
  as_draws_df() |> 
  dplyr::select(matches("^b_.*|^sd_.*")) |> 
  exp() |> 
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat,
    ess_bulk,
    ess_tail,
    Pl = ~mean(.x < 1),
    Pg = ~mean(.x >1)
  )

# Check evidence for Pl and Pg because the limits for power and upper gives 97.5% confidence, Pl and Pg is more fine tuned than that for this case.

# When there is no symbionts there are 13x more likely to be predated on.
# When crabs are present on average odds of predation are 55% less likely. NO evidence on that
# On average shrimp presence declines predation 84%, no sufficient evidence again.
# When both are present, odds of predation decrease by 95% and we have strong evidence on that.

#sd_BLOCK__Intercept -> variability between corals. 4.84 units is quite low. 
```


# Model investigation
```{r}
 mckeon.brm4 |> 
  emmeans(~SYMBIONT, type = "link") |> 
  pairs() |> 
  gather_emmeans_draws() |> 
  mutate(.value = exp(.value)) |> 
  dplyr::select(-.chain) |> 
  summarise_draws(
    median,
    ~HDInterval::hdi(.x),
            Pl = ~mean(.x < 1),
            Pg = ~mean(.x > 1)
            )
```


# Further analyses


# References