Bayesian GLMM Part9

1 Preparations

Load the necessary libraries

library(car) # for regression diagnostics
library(broom) # for tidy output
library(ggfortify) # for model diagnostics
library(knitr) # for kable
library(effects) # for partial effects plots
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(tidyverse) # for data wrangling
library(brms)
library(tidybayes)
library(bayesplot)
library(broom.mixed)
library(rstan)
library(patchwork)
library(modelsummary)
library(DHARMa)
source("helperFunctions.R")

2 Original data

hughes <- read_csv("../data/hughes_full.csv", trim_ws = TRUE)
glimpse(hughes)

Bleaching Score	2016 Equates to	2017 Equates to
0	No bleaching	No bleaching
1	0-10% bleaching	0-10% bleaching
2	10-30% bleaching	10-30% bleaching
3	30-60% bleaching	30-60% bleaching
4	> 60% bleaching	60-80% bleaching
5	-	> 80% bleaching

3 Data processing

• Make a categorical version of Year
• Combine Scores 4 and 5 into a single Score of 4 to be consistent across years
• Make Score an ordered categorical variable
• Generate a numeric version of categorical Score (as required for the ocat family. Note this must be a 1 indexed numeric.
• Relevel Sector such that it is in the order North, Central, South
• Ensure that there is a categorical version of Reef
• Generate a numeric version of the categorical Reef

hughes <- hughes |>
  mutate(
    fYear = factor(Year),
    Score = ifelse(Score == 5, 4, Score),
    oScore = factor(Score, ordered = TRUE),
    nScore = as.numeric(factor(Score, ordered = TRUE)),
    SectorThree = factor(SectorThree, levels = c("North", "Central", "South")),
    fReef = factor(ReefID),
    nReef = as.numeric(fReef)
  )
# now make a version that is just 2016
hughes <- hughes |>
  filter(fYear == 2016) |>
  dplyr::select(REEF = ReefID, HABITAT = Habitat, SECTOR = SectorThree, SCORE = Score)
write_csv(hughes, file = "../data/hughes.csv")
## hughes.colors = c('#FFFFFF', rev(heat.colors(length(levels(hughes$oScore))))[-1])

4 Scenario

Once it is established that mass coral bleaching is occurring on the Great Barrier Reef (GBR), a monitoring team is mobilised in both the air and water in order to document the scale and extent of the bleaching. To better understand the causes and consequences of bleaching, one marine ecologist (lets call him Terry) was interested in investigating differences in bleaching extent across different reef habitats. To do so, aerial surveys were partitioned into four habitats (C - crest, F - flank, L - lower, U - upper).

Bleaching is scored categorically according to the following scale.

Bleaching Score	2016 Equates to
0	No bleaching
1	0-10% bleaching
2	10-30% bleaching
3	30-60% bleaching
4	> 60% bleaching

The GBR is very large and the extent of coral bleaching is not uniform throughout the reef. Hence, Terry wanted to see if the habitat patterns were similar throughout the GBR or whether they were dependent on the overall bleaching severity.

5 Read in the data

2016 data only

hughes <- read_csv("../data/hughes.csv", trim_ws = TRUE)

Rows: 3394 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): REEF, HABITAT, SECTOR
dbl (1): SCORE

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

glimpse(hughes)

Rows: 3,394
Columns: 4
$ REEF    <chr> "09-357", "09-357", "09-357", "09-357", "09-357", "09-366", "0…
$ HABITAT <chr> "C", "C", "F", "F", "U", "C", "L", "C", "C", "L", "C", "C", "U…
$ SECTOR  <chr> "North", "North", "North", "North", "North", "North", "North",…
$ SCORE   <dbl> 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 2,…

REEF	HABITAT	SECTOR	SCORE
09-357	C	North	4
09-357	F	North	4
09-357	U	North	3

glimpse

glimpse(hughes)

Rows: 3,394
Columns: 4
$ REEF    <chr> "09-357", "09-357", "09-357", "09-357", "09-357", "09-366", "0…
$ HABITAT <chr> "C", "C", "F", "F", "U", "C", "L", "C", "C", "L", "C", "C", "U…
$ SECTOR  <chr> "North", "North", "North", "North", "North", "North", "North",…
$ SCORE   <dbl> 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, 2,…

head

## Explore the first 6 rows of the data
head(hughes)

# A tibble: 6 × 4
  REEF   HABITAT SECTOR SCORE
  <chr>  <chr>   <chr>  <dbl>
1 09-357 C       North      4
2 09-357 C       North      4
3 09-357 F       North      4
4 09-357 F       North      4
5 09-357 U       North      4
6 09-366 C       North      3

str

str(hughes)

spc_tbl_ [3,394 × 4] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ REEF   : chr [1:3394] "09-357" "09-357" "09-357" "09-357" ...
 $ HABITAT: chr [1:3394] "C" "C" "F" "F" ...
 $ SECTOR : chr [1:3394] "North" "North" "North" "North" ...
 $ SCORE  : num [1:3394] 4 4 4 4 4 3 3 3 3 3 ...
 - attr(*, "spec")=
  .. cols(
  ..   REEF = col_character(),
  ..   HABITAT = col_character(),
  ..   SECTOR = col_character(),
  ..   SCORE = col_double()
  .. )
 - attr(*, "problems")=<externalptr> 

Easystats (datawizard)

hughes |> datawizard::data_codebook()

hughes (3394 rows and 4 variables, 4 shown)

ID | Name    | Type      | Missings |  Values |            N
---+---------+-----------+----------+---------+-------------
1  | REEF    | character | 0 (0.0%) |  09-304 |    5 ( 0.1%)
   |         |           |          |  09-306 |    2 ( 0.1%)
   |         |           |          |  09-355 |    5 ( 0.1%)
   |         |           |          |  09-357 |    5 ( 0.1%)
   |         |           |          | 09-358a |    3 ( 0.1%)
   |         |           |          |  09-366 |    2 ( 0.1%)
   |         |           |          |  09-368 |    3 ( 0.1%)
   |         |           |          |  09-372 |    2 ( 0.1%)
   |         |           |          |  09-376 |    7 ( 0.2%)
   |         |           |          |  09-378 |    4 ( 0.1%)
   |         |           |          |   (...) |             
---+---------+-----------+----------+---------+-------------
2  | HABITAT | character | 0 (0.0%) |       C | 1919 (56.5%)
   |         |           |          |       F |  646 (19.0%)
   |         |           |          |       L |  604 (17.8%)
   |         |           |          |       U |  225 ( 6.6%)
---+---------+-----------+----------+---------+-------------
3  | SECTOR  | character | 0 (0.0%) | Central |  439 (12.9%)
   |         |           |          |   North | 1546 (45.6%)
   |         |           |          |   South | 1409 (41.5%)
---+---------+-----------+----------+---------+-------------
4  | SCORE   | numeric   | 0 (0.0%) |       0 |  475 (14.0%)
   |         |           |          |       1 |  698 (20.6%)
   |         |           |          |       2 |  641 (18.9%)
   |         |           |          |       3 |  587 (17.3%)
   |         |           |          |       4 |  993 (29.3%)
------------------------------------------------------------

Skim (modelsummary)

hughes |> modelsummary::datasummary_skim(by = "HABITAT")

Warning: These variables were omitted because they include more than 50 levels:
REEF.

	HABITAT	Unique	Missing Pct.	Mean	SD	Min	Median	Max	Histogram
SCORE	C	5	0	2.5	1.4	0.0	3.0	4.0
	F	5	0	2.0	1.4	0.0	2.0	4.0
	U	5	0	2.7	1.3	0.0	3.0	4.0
	L	5	0	1.6	1.4	0.0	1.0	4.0
		N	%
HABITAT	C	1919	56.5
	F	646	19.0
	L	604	17.8
	U	225	6.6
SECTOR	Central	439	12.9
	North	1546	45.6
	South	1409	41.5

6 Data preparation

• ensure score is ordinal
• declare categorical variables as factors - and relevel sector
• declare random effects as factors

hughes <- hughes |>
  mutate(
    oSCORE = factor(SCORE, ordered = TRUE),
    HABITAT = factor(HABITAT),
    SECTOR = factor(SECTOR, levels = c("North", "Central", "South")),
    REEF = factor(REEF)
  )
## hughes.colors = c("#FFFFFF", rev(heat.colors(length(levels(hughes$oSCORE))))[-1])

7 Exploratory Data Analysis

Proportional (cumulative link) odds-ratio models are useful when the latent (unmeasured) response is recorded on a ordinal (ordered categorical) scale. When this is the case, we can calculate the probability of a that an observed ordinal score (\(y\)) is less than or equal to any given level (\(i\): category) given a set of covariates (\(\boldsymbol{X}\)) according to the following:

\[ Pr(y\le i|\boldsymbol{X}) = \frac{1}{1+e^{-(\theta_i - \boldsymbol{\beta}.\boldsymbol{X})}}\\ \]

where \(y\) is the observed categorical response, \(\boldsymbol{X}\) is a (\(n \times p\)) effects model matrix, \(\boldsymbol{\beta}\) are the \(p\) effects parameters and \(\theta_i\) are the \(K-1\) category thresholds

Dot plots

• count the number of samples in each bleaching category per sector per habitat
• express this number as proportions
• reverse the order of the score levels

# Scatterplot
hughes |>
  ggplot(aes(y = oSCORE, x = HABITAT)) +
  geom_point(position = position_jitter()) +
  facet_wrap(~SECTOR)

Line plot

too messy to make much of this…

hughes |>
  group_by(SECTOR, REEF, HABITAT) |>
  summarise(SCORE = mean(SCORE)) |>
  ungroup() |>
  ggplot(aes(y = SCORE, x = as.numeric(HABITAT), group = REEF)) +
  geom_blank(aes(x = HABITAT)) +
  geom_line() +
  facet_grid(~SECTOR)

`summarise()` has grouped output by 'SECTOR', 'REEF'. You can override using
the `.groups` argument.

Stacked bar

hughes |>
  group_by(SECTOR, HABITAT, oSCORE) |>
  summarise(n = n()) |>
  ungroup() |>
  group_by(SECTOR, HABITAT) |>
  mutate(prop = n / sum(n)) |>
  mutate(oSCORE = factor(oSCORE, levels = rev(levels(oSCORE)))) ->
hughes.sum

`summarise()` has grouped output by 'SECTOR', 'HABITAT'. You can override using
the `.groups` argument.

## hughes.sum <- hughes |>
##     count(SECTOR,HABITAT,oSCORE) |>
##     group_by(SECTOR, HABITAT) |>
##     mutate(prop=prop.table(n),
##            oSCORE=factor(oSCORE, levels=rev(levels(oSCORE))))

hughes.sum |> head()

# A tibble: 6 × 5
# Groups:   SECTOR, HABITAT [2]
  SECTOR HABITAT oSCORE     n   prop
  <fct>  <fct>   <ord>  <int>  <dbl>
1 North  C       0         17 0.0172
2 North  C       1         47 0.0475
3 North  C       2        110 0.111 
4 North  C       3        191 0.193 
5 North  C       4        625 0.631 
6 North  F       0          3 0.0121

ggplot(data = hughes.sum, aes(y = prop, x = HABITAT)) +
  geom_bar(stat = "Identity", aes(fill = oSCORE), color = "black") +
  facet_grid(~SECTOR) +
  ## scale_fill_manual('Bleaching score', values=rev(hughes.colors) ) +
  scale_fill_manual("Bleaching score", values = c(heat.colors(5)[-5], "#FFFFFF")) +
  scale_y_continuous("Proportion of Reef", expand = c(0, 0)) +
  theme_bw() +
  theme(panel.spacing.y = unit(10, "pt"))

8 Fit the model

brms

Using default priors

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularisation (to help prevent over fitting) and help stabilise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

hughes.form <- bf(oSCORE ~ HABITAT*SECTOR + (1|REEF),
  family = cumulative(link = "logit", threshold = "flexible")
)
options(width=150)
hughes.form %>% get_prior(data = hughes)

                prior     class                   coef group resp dpar nlpar lb ub       source
               (flat)         b                                                         default
               (flat)         b               HABITATF                             (vectorized)
               (flat)         b HABITATF:SECTORCentral                             (vectorized)
               (flat)         b   HABITATF:SECTORSouth                             (vectorized)
               (flat)         b               HABITATL                             (vectorized)
               (flat)         b HABITATL:SECTORCentral                             (vectorized)
               (flat)         b   HABITATL:SECTORSouth                             (vectorized)
               (flat)         b               HABITATU                             (vectorized)
               (flat)         b HABITATU:SECTORCentral                             (vectorized)
               (flat)         b   HABITATU:SECTORSouth                             (vectorized)
               (flat)         b          SECTORCentral                             (vectorized)
               (flat)         b            SECTORSouth                             (vectorized)
 student_t(3, 0, 2.5) Intercept                                                         default
 student_t(3, 0, 2.5) Intercept                      1                             (vectorized)
 student_t(3, 0, 2.5) Intercept                      2                             (vectorized)
 student_t(3, 0, 2.5) Intercept                      3                             (vectorized)
 student_t(3, 0, 2.5) Intercept                      4                             (vectorized)
 student_t(3, 0, 2.5)        sd                                               0         default
 student_t(3, 0, 2.5)        sd                         REEF                  0    (vectorized)
 student_t(3, 0, 2.5)        sd              Intercept  REEF                  0    (vectorized)

options(width=80)

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucked out of thin air):

• \(\beta_0\): normal centred at 84 with a standard deviation of 5.9
  • mean of 84: since median(starling$MASS)
  • sd of 5.9: since mad(starling$MASS)
• \(\beta_1\): normal centred at 0 with a standard deviation of 13
  • sd of 13: since sd(starling$MASS) / apply(model.matrix(~MONTH*SITUATION, starling)[, -1], 2, sd)
• \(\beta_{2-4}\): normal centred at 0 with a standard deviation of 15
  • sd of 15: since sd(starling$MASS) / apply(model.matrix(~MONTH*SITUATION, starling)[, -1], 2, sd)
• \(\beta_{5-7}\): normal centred at 0 with a standard deviation of 20
  • sd of 20: since sd(starling$MASS) / apply(model.matrix(~MONTH*SITUATION, starling)[, -1], 2, sd)
• \(\sigma\): gamma with shape parameters 6.5 and 1
  • 6.5: since sd(starling$MASS)
• \(\sigma_j\): half-cauchy with parameters 0 and 2.5.
  • 2.: since sqrt(sd(starling$MASS))
• \(\Sigma\): decov with:
  • regularisation: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
  • concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
  • shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.

Note, for hierarchical models, the model will tend to want to have a large \(sigma\) in order to fit the data better. It is a good idea to regularise this tendency by applying a prior that has most mass around zero. Suitable candidates include:

• half-t: as the degrees of freedom approach infinity, this will approach a half-normal
• half-cauchy: this is essentially a half-t with 1 degree of freedom
• exponential

I will also overlay the raw data for comparison.

priors <- prior(normal(0, 1), class = "Intercept") +
  prior(normal(0, 1), class = "b") +
  prior(student_t(3, 0, 1), class = "sd")
hughes.form <- bf(oSCORE ~ HABITAT * SECTOR + (1 | REEF),
  family = cumulative(link = "logit", threshold = "flexible")
)
hughes.brm2 <- brm(hughes.form,
  data = hughes,
  prior = priors,
  sample_prior = "only",
  iter = 5000,
  warmup = 1000,
  chains = 3, cores = 3,
  thin = 5,
  refresh = 0,
  backend = "rstan"
)

Compiling Stan program...

Start sampling

## hughes.brm2 %>%
##     ggpredict(~HABITAT*SECTOR) |>
##     plot(show_data = TRUE)

hughes.brm2 |>
  conditional_effects("HABITAT",
    conditions = make_conditions(hughes.brm2, "SECTOR"),
    categorical = TRUE
  )

The above seem sufficiently wide whilst at the same time not providing any encouragement for the sampler to wander off into very unsupported areas.

hughes.brm3 <- update(hughes.brm2,
  sample_prior = "yes",
  control = list(adapt_delta = 0.99),
  refresh = 100
)

The desired updates require recompiling the model

Compiling Stan program...

Start sampling

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001547 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 15.47 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1: Iteration:  100 / 5000 [  2%]  (Warmup)
Chain 1: Iteration:  200 / 5000 [  4%]  (Warmup)
Chain 1: Iteration:  300 / 5000 [  6%]  (Warmup)
Chain 1: Iteration:  400 / 5000 [  8%]  (Warmup)
Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1: Iteration:  600 / 5000 [ 12%]  (Warmup)
Chain 1: Iteration:  700 / 5000 [ 14%]  (Warmup)
Chain 1: Iteration:  800 / 5000 [ 16%]  (Warmup)
Chain 1: Iteration:  900 / 5000 [ 18%]  (Warmup)
Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1: Iteration: 1001 / 5000 [ 20%]  (Sampling)
Chain 1: Iteration: 1100 / 5000 [ 22%]  (Sampling)
Chain 1: Iteration: 1200 / 5000 [ 24%]  (Sampling)
Chain 1: Iteration: 1300 / 5000 [ 26%]  (Sampling)
Chain 1: Iteration: 1400 / 5000 [ 28%]  (Sampling)
Chain 1: Iteration: 1500 / 5000 [ 30%]  (Sampling)
Chain 1: Iteration: 1600 / 5000 [ 32%]  (Sampling)
Chain 1: Iteration: 1700 / 5000 [ 34%]  (Sampling)
Chain 1: Iteration: 1800 / 5000 [ 36%]  (Sampling)
Chain 1: Iteration: 1900 / 5000 [ 38%]  (Sampling)
Chain 1: Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 1: Iteration: 2100 / 5000 [ 42%]  (Sampling)
Chain 1: Iteration: 2200 / 5000 [ 44%]  (Sampling)
Chain 1: Iteration: 2300 / 5000 [ 46%]  (Sampling)
Chain 1: Iteration: 2400 / 5000 [ 48%]  (Sampling)
Chain 1: Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 1: Iteration: 2600 / 5000 [ 52%]  (Sampling)
Chain 1: Iteration: 2700 / 5000 [ 54%]  (Sampling)
Chain 1: Iteration: 2800 / 5000 [ 56%]  (Sampling)
Chain 1: Iteration: 2900 / 5000 [ 58%]  (Sampling)
Chain 1: Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 1: Iteration: 3100 / 5000 [ 62%]  (Sampling)
Chain 1: Iteration: 3200 / 5000 [ 64%]  (Sampling)
Chain 1: Iteration: 3300 / 5000 [ 66%]  (Sampling)
Chain 1: Iteration: 3400 / 5000 [ 68%]  (Sampling)
Chain 1: Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 1: Iteration: 3600 / 5000 [ 72%]  (Sampling)
Chain 1: Iteration: 3700 / 5000 [ 74%]  (Sampling)
Chain 1: Iteration: 3800 / 5000 [ 76%]  (Sampling)
Chain 1: Iteration: 3900 / 5000 [ 78%]  (Sampling)
Chain 1: Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 1: Iteration: 4100 / 5000 [ 82%]  (Sampling)
Chain 1: Iteration: 4200 / 5000 [ 84%]  (Sampling)
Chain 1: Iteration: 4300 / 5000 [ 86%]  (Sampling)
Chain 1: Iteration: 4400 / 5000 [ 88%]  (Sampling)
Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1: Iteration: 4600 / 5000 [ 92%]  (Sampling)
Chain 1: Iteration: 4700 / 5000 [ 94%]  (Sampling)
Chain 1: Iteration: 4800 / 5000 [ 96%]  (Sampling)
Chain 1: Iteration: 4900 / 5000 [ 98%]  (Sampling)
Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 165.897 seconds (Warm-up)
Chain 1:                759.202 seconds (Sampling)
Chain 1:                925.099 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.001502 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 15.02 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 2: Iteration:  100 / 5000 [  2%]  (Warmup)
Chain 2: Iteration:  200 / 5000 [  4%]  (Warmup)
Chain 2: Iteration:  300 / 5000 [  6%]  (Warmup)
Chain 2: Iteration:  400 / 5000 [  8%]  (Warmup)
Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2: Iteration:  600 / 5000 [ 12%]  (Warmup)
Chain 2: Iteration:  700 / 5000 [ 14%]  (Warmup)
Chain 2: Iteration:  800 / 5000 [ 16%]  (Warmup)
Chain 2: Iteration:  900 / 5000 [ 18%]  (Warmup)
Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2: Iteration: 1001 / 5000 [ 20%]  (Sampling)
Chain 2: Iteration: 1100 / 5000 [ 22%]  (Sampling)
Chain 2: Iteration: 1200 / 5000 [ 24%]  (Sampling)
Chain 2: Iteration: 1300 / 5000 [ 26%]  (Sampling)
Chain 2: Iteration: 1400 / 5000 [ 28%]  (Sampling)
Chain 2: Iteration: 1500 / 5000 [ 30%]  (Sampling)
Chain 2: Iteration: 1600 / 5000 [ 32%]  (Sampling)
Chain 2: Iteration: 1700 / 5000 [ 34%]  (Sampling)
Chain 2: Iteration: 1800 / 5000 [ 36%]  (Sampling)
Chain 2: Iteration: 1900 / 5000 [ 38%]  (Sampling)
Chain 2: Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 2: Iteration: 2100 / 5000 [ 42%]  (Sampling)
Chain 2: Iteration: 2200 / 5000 [ 44%]  (Sampling)
Chain 2: Iteration: 2300 / 5000 [ 46%]  (Sampling)
Chain 2: Iteration: 2400 / 5000 [ 48%]  (Sampling)
Chain 2: Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 2: Iteration: 2600 / 5000 [ 52%]  (Sampling)
Chain 2: Iteration: 2700 / 5000 [ 54%]  (Sampling)
Chain 2: Iteration: 2800 / 5000 [ 56%]  (Sampling)
Chain 2: Iteration: 2900 / 5000 [ 58%]  (Sampling)
Chain 2: Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 2: Iteration: 3100 / 5000 [ 62%]  (Sampling)
Chain 2: Iteration: 3200 / 5000 [ 64%]  (Sampling)
Chain 2: Iteration: 3300 / 5000 [ 66%]  (Sampling)
Chain 2: Iteration: 3400 / 5000 [ 68%]  (Sampling)
Chain 2: Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 2: Iteration: 3600 / 5000 [ 72%]  (Sampling)
Chain 2: Iteration: 3700 / 5000 [ 74%]  (Sampling)
Chain 2: Iteration: 3800 / 5000 [ 76%]  (Sampling)
Chain 2: Iteration: 3900 / 5000 [ 78%]  (Sampling)
Chain 2: Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 2: Iteration: 4100 / 5000 [ 82%]  (Sampling)
Chain 2: Iteration: 4200 / 5000 [ 84%]  (Sampling)
Chain 2: Iteration: 4300 / 5000 [ 86%]  (Sampling)
Chain 2: Iteration: 4400 / 5000 [ 88%]  (Sampling)
Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2: Iteration: 4600 / 5000 [ 92%]  (Sampling)
Chain 2: Iteration: 4700 / 5000 [ 94%]  (Sampling)
Chain 2: Iteration: 4800 / 5000 [ 96%]  (Sampling)
Chain 2: Iteration: 4900 / 5000 [ 98%]  (Sampling)
Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 159.13 seconds (Warm-up)
Chain 2:                531.695 seconds (Sampling)
Chain 2:                690.825 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.001552 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 15.52 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 3: Iteration:  100 / 5000 [  2%]  (Warmup)
Chain 3: Iteration:  200 / 5000 [  4%]  (Warmup)
Chain 3: Iteration:  300 / 5000 [  6%]  (Warmup)
Chain 3: Iteration:  400 / 5000 [  8%]  (Warmup)
Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3: Iteration:  600 / 5000 [ 12%]  (Warmup)
Chain 3: Iteration:  700 / 5000 [ 14%]  (Warmup)
Chain 3: Iteration:  800 / 5000 [ 16%]  (Warmup)
Chain 3: Iteration:  900 / 5000 [ 18%]  (Warmup)
Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3: Iteration: 1001 / 5000 [ 20%]  (Sampling)
Chain 3: Iteration: 1100 / 5000 [ 22%]  (Sampling)
Chain 3: Iteration: 1200 / 5000 [ 24%]  (Sampling)
Chain 3: Iteration: 1300 / 5000 [ 26%]  (Sampling)
Chain 3: Iteration: 1400 / 5000 [ 28%]  (Sampling)
Chain 3: Iteration: 1500 / 5000 [ 30%]  (Sampling)
Chain 3: Iteration: 1600 / 5000 [ 32%]  (Sampling)
Chain 3: Iteration: 1700 / 5000 [ 34%]  (Sampling)
Chain 3: Iteration: 1800 / 5000 [ 36%]  (Sampling)
Chain 3: Iteration: 1900 / 5000 [ 38%]  (Sampling)
Chain 3: Iteration: 2000 / 5000 [ 40%]  (Sampling)
Chain 3: Iteration: 2100 / 5000 [ 42%]  (Sampling)
Chain 3: Iteration: 2200 / 5000 [ 44%]  (Sampling)
Chain 3: Iteration: 2300 / 5000 [ 46%]  (Sampling)
Chain 3: Iteration: 2400 / 5000 [ 48%]  (Sampling)
Chain 3: Iteration: 2500 / 5000 [ 50%]  (Sampling)
Chain 3: Iteration: 2600 / 5000 [ 52%]  (Sampling)
Chain 3: Iteration: 2700 / 5000 [ 54%]  (Sampling)
Chain 3: Iteration: 2800 / 5000 [ 56%]  (Sampling)
Chain 3: Iteration: 2900 / 5000 [ 58%]  (Sampling)
Chain 3: Iteration: 3000 / 5000 [ 60%]  (Sampling)
Chain 3: Iteration: 3100 / 5000 [ 62%]  (Sampling)
Chain 3: Iteration: 3200 / 5000 [ 64%]  (Sampling)
Chain 3: Iteration: 3300 / 5000 [ 66%]  (Sampling)
Chain 3: Iteration: 3400 / 5000 [ 68%]  (Sampling)
Chain 3: Iteration: 3500 / 5000 [ 70%]  (Sampling)
Chain 3: Iteration: 3600 / 5000 [ 72%]  (Sampling)
Chain 3: Iteration: 3700 / 5000 [ 74%]  (Sampling)
Chain 3: Iteration: 3800 / 5000 [ 76%]  (Sampling)
Chain 3: Iteration: 3900 / 5000 [ 78%]  (Sampling)
Chain 3: Iteration: 4000 / 5000 [ 80%]  (Sampling)
Chain 3: Iteration: 4100 / 5000 [ 82%]  (Sampling)
Chain 3: Iteration: 4200 / 5000 [ 84%]  (Sampling)
Chain 3: Iteration: 4300 / 5000 [ 86%]  (Sampling)
Chain 3: Iteration: 4400 / 5000 [ 88%]  (Sampling)
Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3: Iteration: 4600 / 5000 [ 92%]  (Sampling)
Chain 3: Iteration: 4700 / 5000 [ 94%]  (Sampling)
Chain 3: Iteration: 4800 / 5000 [ 96%]  (Sampling)
Chain 3: Iteration: 4900 / 5000 [ 98%]  (Sampling)
Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 181.048 seconds (Warm-up)
Chain 3:                381.868 seconds (Sampling)
Chain 3:                562.916 seconds (Total)
Chain 3: 

##
## load(file = '../ws/testing/hughes.brm3')

hughes.brm3 |>
  conditional_effects("HABITAT", categorical = TRUE) |>
  plot(points = TRUE)

hughes.brm3 |>
  conditional_effects("SECTOR", categorical = TRUE) |>
  plot(points = TRUE)

hughes.brm3 |>
  conditional_effects("HABITAT",
    conditions = make_conditions(hughes.brm3, "SECTOR"),
    categorical = TRUE
  )

a <-
  hughes.brm3 |>
  conditional_effects("HABITAT",
    conditions = make_conditions(hughes.brm3, "SECTOR"),
    categorical = TRUE,
  )

Plotting prior and posterior

hughes.brm3 %>% get_variables()

  [1] "b_Intercept[1]"            "b_Intercept[2]"           
  [3] "b_Intercept[3]"            "b_Intercept[4]"           
  [5] "b_HABITATF"                "b_HABITATL"               
  [7] "b_HABITATU"                "b_SECTORCentral"          
  [9] "b_SECTORSouth"             "b_HABITATF:SECTORCentral" 
 [11] "b_HABITATL:SECTORCentral"  "b_HABITATU:SECTORCentral" 
 [13] "b_HABITATF:SECTORSouth"    "b_HABITATL:SECTORSouth"   
 [15] "b_HABITATU:SECTORSouth"    "sd_REEF__Intercept"       
 [17] "disc"                      "Intercept[1]"             
 [19] "Intercept[2]"              "Intercept[3]"             
 [21] "Intercept[4]"              "r_REEF[09-304,Intercept]" 
 [23] "r_REEF[09-306,Intercept]"  "r_REEF[09-355,Intercept]" 
 [25] "r_REEF[09-357,Intercept]"  "r_REEF[09-358a,Intercept]"
 [27] "r_REEF[09-366,Intercept]"  "r_REEF[09-368,Intercept]" 
 [29] "r_REEF[09-372,Intercept]"  "r_REEF[09-376,Intercept]" 
 [31] "r_REEF[09-378,Intercept]"  "r_REEF[09-380,Intercept]" 
 [33] "r_REEF[09-394,Intercept]"  "r_REEF[09-406,Intercept]" 
 [35] "r_REEF[09-407,Intercept]"  "r_REEF[09-408,Intercept]" 
 [37] "r_REEF[09-409,Intercept]"  "r_REEF[09-410,Intercept]" 
 [39] "r_REEF[09-411,Intercept]"  "r_REEF[09-419a,Intercept]"
 [41] "r_REEF[09-423,Intercept]"  "r_REEF[09-424,Intercept]" 
 [43] "r_REEF[10-283b,Intercept]" "r_REEF[10-299,Intercept]" 
 [45] "r_REEF[10-302,Intercept]"  "r_REEF[10-303,Intercept]" 
 [47] "r_REEF[10-303d,Intercept]" "r_REEF[10-305,Intercept]" 
 [49] "r_REEF[10-318a,Intercept]" "r_REEF[10-318c,Intercept]"
 [51] "r_REEF[10-320,Intercept]"  "r_REEF[10-323,Intercept]" 
 [53] "r_REEF[10-325,Intercept]"  "r_REEF[10-326,Intercept]" 
 [55] "r_REEF[10-327,Intercept]"  "r_REEF[10-335,Intercept]" 
 [57] "r_REEF[10-336,Intercept]"  "r_REEF[10-337,Intercept]" 
 [59] "r_REEF[10-338,Intercept]"  "r_REEF[10-341,Intercept]" 
 [61] "r_REEF[10-347a,Intercept]" "r_REEF[10-349,Intercept]" 
 [63] "r_REEF[10-350,Intercept]"  "r_REEF[10-351,Intercept]" 
 [65] "r_REEF[10-352,Intercept]"  "r_REEF[10-353,Intercept]" 
 [67] "r_REEF[10-355,Intercept]"  "r_REEF[10-377,Intercept]" 
 [69] "r_REEF[10-381,Intercept]"  "r_REEF[10-384,Intercept]" 
 [71] "r_REEF[10-386,Intercept]"  "r_REEF[10-387,Intercept]" 
 [73] "r_REEF[10-390,Intercept]"  "r_REEF[10-391,Intercept]" 
 [75] "r_REEF[10-399,Intercept]"  "r_REEF[10-402,Intercept]" 
 [77] "r_REEF[10-406,Intercept]"  "r_REEF[10-412,Intercept]" 
 [79] "r_REEF[10-414,Intercept]"  "r_REEF[10-416b,Intercept]"
 [81] "r_REEF[10-418a,Intercept]" "r_REEF[10-419a,Intercept]"
 [83] "r_REEF[10-420,Intercept]"  "r_REEF[10-426,Intercept]" 
 [85] "r_REEF[10-441,Intercept]"  "r_REEF[10-449,Intercept]" 
 [87] "r_REEF[10-450,Intercept]"  "r_REEF[10-456a,Intercept]"
 [89] "r_REEF[10-458,Intercept]"  "r_REEF[10-458b,Intercept]"
 [91] "r_REEF[10-503,Intercept]"  "r_REEF[10-514,Intercept]" 
 [93] "r_REEF[10-518c,Intercept]" "r_REEF[10-521a,Intercept]"
 [95] "r_REEF[10-541,Intercept]"  "r_REEF[10-547b,Intercept]"
 [97] "r_REEF[10-562,Intercept]"  "r_REEF[10-565,Intercept]" 
 [99] "r_REEF[10-566a,Intercept]" "r_REEF[10-567,Intercept]" 
[101] "r_REEF[10-572,Intercept]"  "r_REEF[10-575,Intercept]" 
[103] "r_REEF[10-576,Intercept]"  "r_REEF[10-582a,Intercept]"
[105] "r_REEF[10-583,Intercept]"  "r_REEF[10-585,Intercept]" 
[107] "r_REEF[10-598a,Intercept]" "r_REEF[10-619,Intercept]" 
[109] "r_REEF[10-624,Intercept]"  "r_REEF[10-800,Intercept]" 
[111] "r_REEF[11-001,Intercept]"  "r_REEF[11-025,Intercept]" 
[113] "r_REEF[11-026,Intercept]"  "r_REEF[11-035,Intercept]" 
[115] "r_REEF[11-038,Intercept]"  "r_REEF[11-052,Intercept]" 
[117] "r_REEF[11-053,Intercept]"  "r_REEF[11-055,Intercept]" 
[119] "r_REEF[11-061,Intercept]"  "r_REEF[11-062,Intercept]" 
[121] "r_REEF[11-072,Intercept]"  "r_REEF[11-073,Intercept]" 
[123] "r_REEF[11-074,Intercept]"  "r_REEF[11-075,Intercept]" 
[125] "r_REEF[11-076,Intercept]"  "r_REEF[11-077,Intercept]" 
[127] "r_REEF[11-081,Intercept]"  "r_REEF[11-082,Intercept]" 
[129] "r_REEF[11-087b,Intercept]" "r_REEF[11-089b,Intercept]"
[131] "r_REEF[11-091,Intercept]"  "r_REEF[11-094,Intercept]" 
[133] "r_REEF[11-102,Intercept]"  "r_REEF[11-109a,Intercept]"
[135] "r_REEF[11-109b,Intercept]" "r_REEF[11-110,Intercept]" 
[137] "r_REEF[11-113,Intercept]"  "r_REEF[11-114,Intercept]" 
[139] "r_REEF[11-115,Intercept]"  "r_REEF[11-116,Intercept]" 
[141] "r_REEF[11-117,Intercept]"  "r_REEF[11-118a,Intercept]"
[143] "r_REEF[11-119,Intercept]"  "r_REEF[11-120,Intercept]" 
[145] "r_REEF[11-123,Intercept]"  "r_REEF[11-125,Intercept]" 
[147] "r_REEF[11-127b,Intercept]" "r_REEF[11-128,Intercept]" 
[149] "r_REEF[11-129,Intercept]"  "r_REEF[11-134,Intercept]" 
[151] "r_REEF[11-136,Intercept]"  "r_REEF[11-137,Intercept]" 
[153] "r_REEF[11-138,Intercept]"  "r_REEF[11-139,Intercept]" 
[155] "r_REEF[11-140,Intercept]"  "r_REEF[11-141,Intercept]" 
[157] "r_REEF[11-142,Intercept]"  "r_REEF[11-143,Intercept]" 
[159] "r_REEF[11-147,Intercept]"  "r_REEF[11-150,Intercept]" 
[161] "r_REEF[11-151,Intercept]"  "r_REEF[11-156,Intercept]" 
[163] "r_REEF[11-161,Intercept]"  "r_REEF[11-173,Intercept]" 
[165] "r_REEF[11-181,Intercept]"  "r_REEF[11-184a,Intercept]"
[167] "r_REEF[11-184b,Intercept]" "r_REEF[11-184c,Intercept]"
[169] "r_REEF[11-187,Intercept]"  "r_REEF[11-188,Intercept]" 
[171] "r_REEF[11-189,Intercept]"  "r_REEF[11-190a,Intercept]"
[173] "r_REEF[11-190b,Intercept]" "r_REEF[11-191,Intercept]" 
[175] "r_REEF[11-194,Intercept]"  "r_REEF[11-197,Intercept]" 
[177] "r_REEF[11-201,Intercept]"  "r_REEF[11-202,Intercept]" 
[179] "r_REEF[11-204,Intercept]"  "r_REEF[11-208a,Intercept]"
[181] "r_REEF[11-210a,Intercept]" "r_REEF[11-212,Intercept]" 
[183] "r_REEF[11-222a,Intercept]" "r_REEF[11-223,Intercept]" 
[185] "r_REEF[11-225,Intercept]"  "r_REEF[11-228,Intercept]" 
[187] "r_REEF[11-229a,Intercept]" "r_REEF[11-229b,Intercept]"
[189] "r_REEF[11-231,Intercept]"  "r_REEF[11-232,Intercept]" 
[191] "r_REEF[11-236,Intercept]"  "r_REEF[11-237a,Intercept]"
[193] "r_REEF[11-238,Intercept]"  "r_REEF[11-239,Intercept]" 
[195] "r_REEF[11-241,Intercept]"  "r_REEF[11-244c,Intercept]"
[197] "r_REEF[12-001,Intercept]"  "r_REEF[12-003a,Intercept]"
[199] "r_REEF[12-004,Intercept]"  "r_REEF[12-006,Intercept]" 
[201] "r_REEF[12-007,Intercept]"  "r_REEF[12-010,Intercept]" 
[203] "r_REEF[12-011,Intercept]"  "r_REEF[12-012,Intercept]" 
[205] "r_REEF[12-013,Intercept]"  "r_REEF[12-014,Intercept]" 
[207] "r_REEF[12-016a,Intercept]" "r_REEF[12-016b,Intercept]"
[209] "r_REEF[12-018,Intercept]"  "r_REEF[12-022,Intercept]" 
[211] "r_REEF[12-027a,Intercept]" "r_REEF[12-031,Intercept]" 
[213] "r_REEF[12-032,Intercept]"  "r_REEF[12-033,Intercept]" 
[215] "r_REEF[12-034,Intercept]"  "r_REEF[12-035,Intercept]" 
[217] "r_REEF[12-036,Intercept]"  "r_REEF[12-037,Intercept]" 
[219] "r_REEF[12-038a,Intercept]" "r_REEF[12-038c,Intercept]"
[221] "r_REEF[12-039,Intercept]"  "r_REEF[12-040,Intercept]" 
[223] "r_REEF[12-042,Intercept]"  "r_REEF[12-043a,Intercept]"
[225] "r_REEF[12-043b,Intercept]" "r_REEF[12-046,Intercept]" 
[227] "r_REEF[12-048,Intercept]"  "r_REEF[12-051,Intercept]" 
[229] "r_REEF[12-053,Intercept]"  "r_REEF[12-056,Intercept]" 
[231] "r_REEF[12-058,Intercept]"  "r_REEF[12-059,Intercept]" 
[233] "r_REEF[12-061,Intercept]"  "r_REEF[12-068,Intercept]" 
[235] "r_REEF[12-069,Intercept]"  "r_REEF[12-070,Intercept]" 
[237] "r_REEF[12-071,Intercept]"  "r_REEF[12-075,Intercept]" 
[239] "r_REEF[12-076,Intercept]"  "r_REEF[12-077,Intercept]" 
[241] "r_REEF[12-078,Intercept]"  "r_REEF[12-094,Intercept]" 
[243] "r_REEF[12-096,Intercept]"  "r_REEF[12-098,Intercept]" 
[245] "r_REEF[12-100,Intercept]"  "r_REEF[12-101,Intercept]" 
[247] "r_REEF[12-102,Intercept]"  "r_REEF[12-104,Intercept]" 
[249] "r_REEF[12-105,Intercept]"  "r_REEF[12-107,Intercept]" 
[251] "r_REEF[12-110,Intercept]"  "r_REEF[12-111,Intercept]" 
[253] "r_REEF[12-112,Intercept]"  "r_REEF[12-114,Intercept]" 
[255] "r_REEF[12-115,Intercept]"  "r_REEF[12-116,Intercept]" 
[257] "r_REEF[12-118,Intercept]"  "r_REEF[12-119,Intercept]" 
[259] "r_REEF[12-127,Intercept]"  "r_REEF[12-129,Intercept]" 
[261] "r_REEF[12-131,Intercept]"  "r_REEF[12-134,Intercept]" 
[263] "r_REEF[12-137,Intercept]"  "r_REEF[12-143,Intercept]" 
[265] "r_REEF[12-145,Intercept]"  "r_REEF[13-002,Intercept]" 
[267] "r_REEF[13-005,Intercept]"  "r_REEF[13-006,Intercept]" 
[269] "r_REEF[13-007,Intercept]"  "r_REEF[13-012,Intercept]" 
[271] "r_REEF[13-014,Intercept]"  "r_REEF[13-019,Intercept]" 
[273] "r_REEF[13-022,Intercept]"  "r_REEF[13-028,Intercept]" 
[275] "r_REEF[13-029,Intercept]"  "r_REEF[13-031,Intercept]" 
[277] "r_REEF[13-032,Intercept]"  "r_REEF[13-034,Intercept]" 
[279] "r_REEF[13-040,Intercept]"  "r_REEF[13-041,Intercept]" 
[281] "r_REEF[13-045,Intercept]"  "r_REEF[13-050a,Intercept]"
[283] "r_REEF[13-055,Intercept]"  "r_REEF[13-056,Intercept]" 
[285] "r_REEF[13-060,Intercept]"  "r_REEF[13-061a,Intercept]"
[287] "r_REEF[13-061b,Intercept]" "r_REEF[13-063,Intercept]" 
[289] "r_REEF[13-069,Intercept]"  "r_REEF[13-072,Intercept]" 
[291] "r_REEF[13-074,Intercept]"  "r_REEF[13-075,Intercept]" 
[293] "r_REEF[13-076,Intercept]"  "r_REEF[13-077,Intercept]" 
[295] "r_REEF[13-078,Intercept]"  "r_REEF[13-079c,Intercept]"
[297] "r_REEF[13-081,Intercept]"  "r_REEF[13-083,Intercept]" 
[299] "r_REEF[13-087,Intercept]"  "r_REEF[13-088,Intercept]" 
[301] "r_REEF[13-091,Intercept]"  "r_REEF[13-093b,Intercept]"
[303] "r_REEF[13-093c,Intercept]" "r_REEF[13-097,Intercept]" 
[305] "r_REEF[13-108,Intercept]"  "r_REEF[13-111,Intercept]" 
[307] "r_REEF[13-116,Intercept]"  "r_REEF[13-117b,Intercept]"
[309] "r_REEF[13-118a,Intercept]" "r_REEF[13-119,Intercept]" 
[311] "r_REEF[13-120,Intercept]"  "r_REEF[13-121,Intercept]" 
[313] "r_REEF[13-122,Intercept]"  "r_REEF[13-125,Intercept]" 
[315] "r_REEF[13-127,Intercept]"  "r_REEF[13-129,Intercept]" 
[317] "r_REEF[13-130,Intercept]"  "r_REEF[13-131,Intercept]" 
[319] "r_REEF[13-133,Intercept]"  "r_REEF[14-003,Intercept]" 
[321] "r_REEF[14-006,Intercept]"  "r_REEF[14-008,Intercept]" 
[323] "r_REEF[14-012,Intercept]"  "r_REEF[14-013,Intercept]" 
[325] "r_REEF[14-016,Intercept]"  "r_REEF[14-017,Intercept]" 
[327] "r_REEF[14-022,Intercept]"  "r_REEF[14-026a,Intercept]"
[329] "r_REEF[14-026b,Intercept]" "r_REEF[14-028,Intercept]" 
[331] "r_REEF[14-029,Intercept]"  "r_REEF[14-032,Intercept]" 
[333] "r_REEF[14-033a,Intercept]" "r_REEF[14-034,Intercept]" 
[335] "r_REEF[14-038,Intercept]"  "r_REEF[14-039,Intercept]" 
[337] "r_REEF[14-042a,Intercept]" "r_REEF[14-045,Intercept]" 
[339] "r_REEF[14-047,Intercept]"  "r_REEF[14-049,Intercept]" 
[341] "r_REEF[14-051,Intercept]"  "r_REEF[14-052,Intercept]" 
[343] "r_REEF[14-053,Intercept]"  "r_REEF[14-055,Intercept]" 
[345] "r_REEF[14-061,Intercept]"  "r_REEF[14-063,Intercept]" 
[347] "r_REEF[14-064,Intercept]"  "r_REEF[14-065,Intercept]" 
[349] "r_REEF[14-066,Intercept]"  "r_REEF[14-068,Intercept]" 
[351] "r_REEF[14-070,Intercept]"  "r_REEF[14-071a,Intercept]"
[353] "r_REEF[14-072,Intercept]"  "r_REEF[14-073,Intercept]" 
[355] "r_REEF[14-074,Intercept]"  "r_REEF[14-075,Intercept]" 
[357] "r_REEF[14-077a,Intercept]" "r_REEF[14-078,Intercept]" 
[359] "r_REEF[14-079,Intercept]"  "r_REEF[14-082,Intercept]" 
[361] "r_REEF[14-085,Intercept]"  "r_REEF[14-086,Intercept]" 
[363] "r_REEF[14-088,Intercept]"  "r_REEF[14-089,Intercept]" 
[365] "r_REEF[14-090,Intercept]"  "r_REEF[14-091,Intercept]" 
[367] "r_REEF[14-093,Intercept]"  "r_REEF[14-094,Intercept]" 
[369] "r_REEF[14-096,Intercept]"  "r_REEF[14-097,Intercept]" 
[371] "r_REEF[14-099,Intercept]"  "r_REEF[14-101,Intercept]" 
[373] "r_REEF[14-102,Intercept]"  "r_REEF[14-103,Intercept]" 
[375] "r_REEF[14-104,Intercept]"  "r_REEF[14-106,Intercept]" 
[377] "r_REEF[14-116a,Intercept]" "r_REEF[14-116b,Intercept]"
[379] "r_REEF[14-116c,Intercept]" "r_REEF[14-116d,Intercept]"
[381] "r_REEF[14-120d,Intercept]" "r_REEF[14-122a,Intercept]"
[383] "r_REEF[14-123,Intercept]"  "r_REEF[14-126,Intercept]" 
[385] "r_REEF[14-132a,Intercept]" "r_REEF[14-132b,Intercept]"
[387] "r_REEF[14-133,Intercept]"  "r_REEF[14-135,Intercept]" 
[389] "r_REEF[14-137,Intercept]"  "r_REEF[14-138,Intercept]" 
[391] "r_REEF[14-139,Intercept]"  "r_REEF[14-140,Intercept]" 
[393] "r_REEF[14-143,Intercept]"  "r_REEF[14-146,Intercept]" 
[395] "r_REEF[14-147a,Intercept]" "r_REEF[14-147b,Intercept]"
[397] "r_REEF[14-150,Intercept]"  "r_REEF[14-154,Intercept]" 
[399] "r_REEF[15-002,Intercept]"  "r_REEF[15-005,Intercept]" 
[401] "r_REEF[15-006,Intercept]"  "r_REEF[15-008,Intercept]" 
[403] "r_REEF[15-009,Intercept]"  "r_REEF[15-012,Intercept]" 
[405] "r_REEF[15-013,Intercept]"  "r_REEF[15-017a,Intercept]"
[407] "r_REEF[15-018,Intercept]"  "r_REEF[15-019,Intercept]" 
[409] "r_REEF[15-021a,Intercept]" "r_REEF[15-022,Intercept]" 
[411] "r_REEF[15-024,Intercept]"  "r_REEF[15-026,Intercept]" 
[413] "r_REEF[15-028,Intercept]"  "r_REEF[15-032,Intercept]" 
[415] "r_REEF[15-033,Intercept]"  "r_REEF[15-034,Intercept]" 
[417] "r_REEF[15-038,Intercept]"  "r_REEF[15-039,Intercept]" 
[419] "r_REEF[15-043,Intercept]"  "r_REEF[15-046,Intercept]" 
[421] "r_REEF[15-049,Intercept]"  "r_REEF[15-050,Intercept]" 
[423] "r_REEF[15-052,Intercept]"  "r_REEF[15-063,Intercept]" 
[425] "r_REEF[15-064,Intercept]"  "r_REEF[15-065,Intercept]" 
[427] "r_REEF[15-066,Intercept]"  "r_REEF[15-071a,Intercept]"
[429] "r_REEF[15-075a,Intercept]" "r_REEF[15-078,Intercept]" 
[431] "r_REEF[15-080,Intercept]"  "r_REEF[15-085,Intercept]" 
[433] "r_REEF[15-086,Intercept]"  "r_REEF[15-087,Intercept]" 
[435] "r_REEF[15-088,Intercept]"  "r_REEF[15-089,Intercept]" 
[437] "r_REEF[15-090,Intercept]"  "r_REEF[15-092,Intercept]" 
[439] "r_REEF[15-093,Intercept]"  "r_REEF[15-094,Intercept]" 
[441] "r_REEF[15-095,Intercept]"  "r_REEF[15-096,Intercept]" 
[443] "r_REEF[15-098,Intercept]"  "r_REEF[15-099a,Intercept]"
[445] "r_REEF[15-099d,Intercept]" "r_REEF[16-015,Intercept]" 
[447] "r_REEF[16-019,Intercept]"  "r_REEF[16-020a,Intercept]"
[449] "r_REEF[16-023,Intercept]"  "r_REEF[16-025,Intercept]" 
[451] "r_REEF[16-026,Intercept]"  "r_REEF[16-029,Intercept]" 
[453] "r_REEF[16-030,Intercept]"  "r_REEF[16-032,Intercept]" 
[455] "r_REEF[16-040,Intercept]"  "r_REEF[16-043a,Intercept]"
[457] "r_REEF[16-046,Intercept]"  "r_REEF[16-049,Intercept]" 
[459] "r_REEF[16-054b,Intercept]" "r_REEF[16-054c,Intercept]"
[461] "r_REEF[16-057,Intercept]"  "r_REEF[16-060,Intercept]" 
[463] "r_REEF[16-068,Intercept]"  "r_REEF[16-071,Intercept]" 
[465] "r_REEF[16-073,Intercept]"  "r_REEF[17-001a,Intercept]"
[467] "r_REEF[17-004,Intercept]"  "r_REEF[17-005,Intercept]" 
[469] "r_REEF[17-006,Intercept]"  "r_REEF[17-008,Intercept]" 
[471] "r_REEF[17-010,Intercept]"  "r_REEF[17-011,Intercept]" 
[473] "r_REEF[17-014,Intercept]"  "r_REEF[17-016a,Intercept]"
[475] "r_REEF[17-017,Intercept]"  "r_REEF[17-018a,Intercept]"
[477] "r_REEF[17-023,Intercept]"  "r_REEF[17-024,Intercept]" 
[479] "r_REEF[17-034,Intercept]"  "r_REEF[17-035,Intercept]" 
[481] "r_REEF[17-037,Intercept]"  "r_REEF[17-042,Intercept]" 
[483] "r_REEF[17-044,Intercept]"  "r_REEF[17-047,Intercept]" 
[485] "r_REEF[17-051,Intercept]"  "r_REEF[17-057,Intercept]" 
[487] "r_REEF[17-059a,Intercept]" "r_REEF[17-062,Intercept]" 
[489] "r_REEF[17-063a,Intercept]" "r_REEF[17-064,Intercept]" 
[491] "r_REEF[17-065,Intercept]"  "r_REEF[17-066,Intercept]" 
[493] "r_REEF[17-067,Intercept]"  "r_REEF[17-068,Intercept]" 
[495] "r_REEF[17-069,Intercept]"  "r_REEF[18-017,Intercept]" 
[497] "r_REEF[18-018,Intercept]"  "r_REEF[18-019,Intercept]" 
[499] "r_REEF[18-023,Intercept]"  "r_REEF[18-024,Intercept]" 
[501] "r_REEF[18-026,Intercept]"  "r_REEF[18-027,Intercept]" 
[503] "r_REEF[18-029,Intercept]"  "r_REEF[18-030,Intercept]" 
[505] "r_REEF[18-031,Intercept]"  "r_REEF[18-032,Intercept]" 
[507] "r_REEF[18-033,Intercept]"  "r_REEF[18-034,Intercept]" 
[509] "r_REEF[18-037,Intercept]"  "r_REEF[18-039,Intercept]" 
[511] "r_REEF[18-042,Intercept]"  "r_REEF[18-043,Intercept]" 
[513] "r_REEF[18-048,Intercept]"  "r_REEF[18-049a,Intercept]"
[515] "r_REEF[18-049b,Intercept]" "r_REEF[18-049c,Intercept]"
[517] "r_REEF[18-049e,Intercept]" "r_REEF[18-054b,Intercept]"
[519] "r_REEF[18-054d,Intercept]" "r_REEF[18-054e,Intercept]"
[521] "r_REEF[18-054f,Intercept]" "r_REEF[18-065,Intercept]" 
[523] "r_REEF[18-069a,Intercept]" "r_REEF[18-069b,Intercept]"
[525] "r_REEF[18-070,Intercept]"  "r_REEF[18-071,Intercept]" 
[527] "r_REEF[18-075,Intercept]"  "r_REEF[18-077,Intercept]" 
[529] "r_REEF[18-078,Intercept]"  "r_REEF[18-084a,Intercept]"
[531] "r_REEF[18-084b,Intercept]" "r_REEF[18-085a,Intercept]"
[533] "r_REEF[18-086,Intercept]"  "r_REEF[18-091,Intercept]" 
[535] "r_REEF[18-096,Intercept]"  "r_REEF[18-099,Intercept]" 
[537] "r_REEF[18-100a,Intercept]" "r_REEF[18-100b,Intercept]"
[539] "r_REEF[18-101,Intercept]"  "r_REEF[18-103,Intercept]" 
[541] "r_REEF[18-105a,Intercept]" "r_REEF[18-105b,Intercept]"
[543] "r_REEF[18-106,Intercept]"  "r_REEF[18-115,Intercept]" 
[545] "r_REEF[18-118,Intercept]"  "r_REEF[18-120,Intercept]" 
[547] "r_REEF[19-009a,Intercept]" "r_REEF[19-009g,Intercept]"
[549] "r_REEF[19-009j,Intercept]" "r_REEF[19-019,Intercept]" 
[551] "r_REEF[19-024,Intercept]"  "r_REEF[19-038a,Intercept]"
[553] "r_REEF[19-038b,Intercept]" "r_REEF[19-038f,Intercept]"
[555] "r_REEF[19-043,Intercept]"  "r_REEF[19-045,Intercept]" 
[557] "r_REEF[19-048,Intercept]"  "r_REEF[19-049,Intercept]" 
[559] "r_REEF[19-050,Intercept]"  "r_REEF[19-051,Intercept]" 
[561] "r_REEF[19-054,Intercept]"  "r_REEF[19-059,Intercept]" 
[563] "r_REEF[19-061,Intercept]"  "r_REEF[19-063a,Intercept]"
[565] "r_REEF[19-065,Intercept]"  "r_REEF[19-066,Intercept]" 
[567] "r_REEF[19-067,Intercept]"  "r_REEF[19-071,Intercept]" 
[569] "r_REEF[19-072a,Intercept]" "r_REEF[19-072b,Intercept]"
[571] "r_REEF[19-072e,Intercept]" "r_REEF[19-074b,Intercept]"
[573] "r_REEF[19-076,Intercept]"  "r_REEF[19-081,Intercept]" 
[575] "r_REEF[19-082,Intercept]"  "r_REEF[19-091,Intercept]" 
[577] "r_REEF[19-096,Intercept]"  "r_REEF[19-097,Intercept]" 
[579] "r_REEF[19-098,Intercept]"  "r_REEF[19-099,Intercept]" 
[581] "r_REEF[19-107,Intercept]"  "r_REEF[19-109a,Intercept]"
[583] "r_REEF[19-113,Intercept]"  "r_REEF[19-114,Intercept]" 
[585] "r_REEF[19-116,Intercept]"  "r_REEF[19-123,Intercept]" 
[587] "r_REEF[19-127,Intercept]"  "r_REEF[19-128,Intercept]" 
[589] "r_REEF[19-129,Intercept]"  "r_REEF[19-134,Intercept]" 
[591] "r_REEF[19-135,Intercept]"  "r_REEF[19-136a,Intercept]"
[593] "r_REEF[19-139,Intercept]"  "r_REEF[19-141,Intercept]" 
[595] "r_REEF[19-146,Intercept]"  "r_REEF[19-147,Intercept]" 
[597] "r_REEF[19-151,Intercept]"  "r_REEF[19-152,Intercept]" 
[599] "r_REEF[19-157,Intercept]"  "r_REEF[19-159,Intercept]" 
[601] "r_REEF[19-167,Intercept]"  "r_REEF[19-168,Intercept]" 
[603] "r_REEF[19-169,Intercept]"  "r_REEF[19-187,Intercept]" 
[605] "r_REEF[19-197,Intercept]"  "r_REEF[19-204,Intercept]" 
[607] "r_REEF[19-206,Intercept]"  "r_REEF[19-207,Intercept]" 
[609] "r_REEF[19-208,Intercept]"  "r_REEF[19-209,Intercept]" 
[611] "r_REEF[19-212,Intercept]"  "r_REEF[19-214,Intercept]" 
[613] "r_REEF[19-216,Intercept]"  "r_REEF[19-220,Intercept]" 
[615] "r_REEF[19-221,Intercept]"  "r_REEF[19-222,Intercept]" 
[617] "r_REEF[19-225,Intercept]"  "r_REEF[19-227,Intercept]" 
[619] "r_REEF[19-246,Intercept]"  "r_REEF[20-014,Intercept]" 
[621] "r_REEF[20-017,Intercept]"  "r_REEF[20-028h,Intercept]"
[623] "r_REEF[20-078a,Intercept]" "r_REEF[20-104,Intercept]" 
[625] "r_REEF[20-105,Intercept]"  "r_REEF[20-112,Intercept]" 
[627] "r_REEF[20-115,Intercept]"  "r_REEF[20-119,Intercept]" 
[629] "r_REEF[20-121,Intercept]"  "r_REEF[20-125,Intercept]" 
[631] "r_REEF[20-126,Intercept]"  "r_REEF[20-127,Intercept]" 
[633] "r_REEF[20-128,Intercept]"  "r_REEF[20-129,Intercept]" 
[635] "r_REEF[20-131,Intercept]"  "r_REEF[20-132,Intercept]" 
[637] "r_REEF[20-134,Intercept]"  "r_REEF[20-135,Intercept]" 
[639] "r_REEF[20-136,Intercept]"  "r_REEF[20-137,Intercept]" 
[641] "r_REEF[20-138,Intercept]"  "r_REEF[20-167,Intercept]" 
[643] "r_REEF[20-172,Intercept]"  "r_REEF[20-179,Intercept]" 
[645] "r_REEF[20-182,Intercept]"  "r_REEF[20-185a,Intercept]"
[647] "r_REEF[20-185c,Intercept]" "r_REEF[20-186,Intercept]" 
[649] "r_REEF[20-187,Intercept]"  "r_REEF[20-194a,Intercept]"
[651] "r_REEF[20-197,Intercept]"  "r_REEF[20-269,Intercept]" 
[653] "r_REEF[20-291,Intercept]"  "r_REEF[20-292,Intercept]" 
[655] "r_REEF[20-295,Intercept]"  "r_REEF[20-296,Intercept]" 
[657] "r_REEF[20-298,Intercept]"  "r_REEF[20-299a,Intercept]"
[659] "r_REEF[20-299c,Intercept]" "r_REEF[20-301a,Intercept]"
[661] "r_REEF[20-301b,Intercept]" "r_REEF[20-302,Intercept]" 
[663] "r_REEF[20-304,Intercept]"  "r_REEF[20-305,Intercept]" 
[665] "r_REEF[20-306,Intercept]"  "r_REEF[20-307,Intercept]" 
[667] "r_REEF[20-308,Intercept]"  "r_REEF[20-322,Intercept]" 
[669] "r_REEF[20-324,Intercept]"  "r_REEF[20-326,Intercept]" 
[671] "r_REEF[20-335,Intercept]"  "r_REEF[20-339a,Intercept]"
[673] "r_REEF[20-339b,Intercept]" "r_REEF[20-348,Intercept]" 
[675] "r_REEF[20-353,Intercept]"  "r_REEF[20-369a,Intercept]"
[677] "r_REEF[20-370a,Intercept]" "r_REEF[20-370b,Intercept]"
[679] "r_REEF[20-377,Intercept]"  "r_REEF[20-378a,Intercept]"
[681] "r_REEF[20-380,Intercept]"  "r_REEF[20-382,Intercept]" 
[683] "r_REEF[20-383,Intercept]"  "r_REEF[20-385,Intercept]" 
[685] "r_REEF[20-386,Intercept]"  "r_REEF[20-388a,Intercept]"
[687] "r_REEF[20-389,Intercept]"  "r_REEF[20-391,Intercept]" 
[689] "r_REEF[20-469,Intercept]"  "r_REEF[20-470,Intercept]" 
[691] "r_REEF[20-474,Intercept]"  "r_REEF[20-568,Intercept]" 
[693] "r_REEF[21-025,Intercept]"  "r_REEF[21-050,Intercept]" 
[695] "r_REEF[21-054,Intercept]"  "r_REEF[21-056b,Intercept]"
[697] "r_REEF[21-057,Intercept]"  "r_REEF[21-061,Intercept]" 
[699] "r_REEF[21-062,Intercept]"  "r_REEF[21-064,Intercept]" 
[701] "r_REEF[21-065,Intercept]"  "r_REEF[21-068,Intercept]" 
[703] "r_REEF[21-070,Intercept]"  "r_REEF[21-071,Intercept]" 
[705] "r_REEF[21-072,Intercept]"  "r_REEF[21-074,Intercept]" 
[707] "r_REEF[21-076,Intercept]"  "r_REEF[21-081,Intercept]" 
[709] "r_REEF[21-082,Intercept]"  "r_REEF[21-086,Intercept]" 
[711] "r_REEF[21-087,Intercept]"  "r_REEF[21-091,Intercept]" 
[713] "r_REEF[21-092,Intercept]"  "r_REEF[21-093,Intercept]" 
[715] "r_REEF[21-094,Intercept]"  "r_REEF[21-100,Intercept]" 
[717] "r_REEF[21-101,Intercept]"  "r_REEF[21-109,Intercept]" 
[719] "r_REEF[21-115,Intercept]"  "r_REEF[21-116,Intercept]" 
[721] "r_REEF[21-128,Intercept]"  "r_REEF[21-129,Intercept]" 
[723] "r_REEF[21-131,Intercept]"  "r_REEF[21-132,Intercept]" 
[725] "r_REEF[21-133,Intercept]"  "r_REEF[21-134,Intercept]" 
[727] "r_REEF[21-135,Intercept]"  "r_REEF[21-136,Intercept]" 
[729] "r_REEF[21-139,Intercept]"  "r_REEF[21-143,Intercept]" 
[731] "r_REEF[21-144,Intercept]"  "r_REEF[21-145a,Intercept]"
[733] "r_REEF[21-146,Intercept]"  "r_REEF[21-147,Intercept]" 
[735] "r_REEF[21-148,Intercept]"  "r_REEF[21-149,Intercept]" 
[737] "r_REEF[21-150,Intercept]"  "r_REEF[21-151,Intercept]" 
[739] "r_REEF[21-152,Intercept]"  "r_REEF[21-157,Intercept]" 
[741] "r_REEF[21-158,Intercept]"  "r_REEF[21-159,Intercept]" 
[743] "r_REEF[21-161,Intercept]"  "r_REEF[21-165,Intercept]" 
[745] "r_REEF[21-172,Intercept]"  "r_REEF[21-180,Intercept]" 
[747] "r_REEF[21-181,Intercept]"  "r_REEF[21-185,Intercept]" 
[749] "r_REEF[21-236,Intercept]"  "r_REEF[21-237,Intercept]" 
[751] "r_REEF[21-238,Intercept]"  "r_REEF[21-241,Intercept]" 
[753] "r_REEF[21-242a,Intercept]" "r_REEF[21-247,Intercept]" 
[755] "r_REEF[21-258,Intercept]"  "r_REEF[21-263,Intercept]" 
[757] "r_REEF[21-289a,Intercept]" "r_REEF[21-290,Intercept]" 
[759] "r_REEF[21-294,Intercept]"  "r_REEF[21-295,Intercept]" 
[761] "r_REEF[21-296,Intercept]"  "r_REEF[21-299,Intercept]" 
[763] "r_REEF[21-300a,Intercept]" "r_REEF[21-303,Intercept]" 
[765] "r_REEF[21-304,Intercept]"  "r_REEF[21-306,Intercept]" 
[767] "r_REEF[21-334,Intercept]"  "r_REEF[21-386,Intercept]" 
[769] "r_REEF[21-394,Intercept]"  "r_REEF[21-403b,Intercept]"
[771] "r_REEF[21-433,Intercept]"  "r_REEF[21-449,Intercept]" 
[773] "r_REEF[21-450,Intercept]"  "r_REEF[21-453,Intercept]" 
[775] "r_REEF[21-460,Intercept]"  "r_REEF[21-464,Intercept]" 
[777] "r_REEF[21-465,Intercept]"  "r_REEF[21-466,Intercept]" 
[779] "r_REEF[21-467,Intercept]"  "r_REEF[21-470,Intercept]" 
[781] "r_REEF[21-474,Intercept]"  "r_REEF[21-479,Intercept]" 
[783] "r_REEF[21-491,Intercept]"  "r_REEF[21-496,Intercept]" 
[785] "r_REEF[21-526,Intercept]"  "r_REEF[21-529,Intercept]" 
[787] "r_REEF[21-530,Intercept]"  "r_REEF[21-531,Intercept]" 
[789] "r_REEF[21-532,Intercept]"  "r_REEF[21-538,Intercept]" 
[791] "r_REEF[21-540,Intercept]"  "r_REEF[21-544,Intercept]" 
[793] "r_REEF[21-546,Intercept]"  "r_REEF[21-551,Intercept]" 
[795] "r_REEF[21-556,Intercept]"  "r_REEF[21-560,Intercept]" 
[797] "r_REEF[21-562,Intercept]"  "r_REEF[21-565,Intercept]" 
[799] "r_REEF[21-566,Intercept]"  "r_REEF[21-570,Intercept]" 
[801] "r_REEF[21-575,Intercept]"  "r_REEF[21-579,Intercept]" 
[803] "r_REEF[21-582,Intercept]"  "r_REEF[21-586,Intercept]" 
[805] "r_REEF[21-591,Intercept]"  "r_REEF[21-761,Intercept]" 
[807] "r_REEF[22-084,Intercept]"  "r_REEF[22-086,Intercept]" 
[809] "r_REEF[22-090,Intercept]"  "r_REEF[22-092,Intercept]" 
[811] "r_REEF[22-093,Intercept]"  "r_REEF[22-103,Intercept]" 
[813] "r_REEF[22-106,Intercept]"  "r_REEF[22-115,Intercept]" 
[815] "r_REEF[22-116,Intercept]"  "r_REEF[22-121,Intercept]" 
[817] "r_REEF[22-122,Intercept]"  "r_REEF[22-123,Intercept]" 
[819] "r_REEF[22-124,Intercept]"  "r_REEF[22-125,Intercept]" 
[821] "r_REEF[22-129,Intercept]"  "r_REEF[22-131,Intercept]" 
[823] "r_REEF[22-133,Intercept]"  "r_REEF[22-134,Intercept]" 
[825] "r_REEF[22-139,Intercept]"  "r_REEF[22-140,Intercept]" 
[827] "r_REEF[22-141,Intercept]"  "r_REEF[22-142,Intercept]" 
[829] "r_REEF[22-143,Intercept]"  "r_REEF[22-144,Intercept]" 
[831] "prior_Intercept"           "prior_b"                  
[833] "prior_sd_REEF"             "lprior"                   
[835] "lp__"                      "accept_stat__"            
[837] "stepsize__"                "treedepth__"              
[839] "n_leapfrog__"              "divergent__"              
[841] "energy__"                 

hughes.brm3 %>%
  hypothesis("HABITATF=0") %>%
  plot()

hughes.brm3 %>%
  hypothesis("SECTORCentral=0") %>%
  plot()

# hughes.brm3 %>% SUYR_prior_and_posterior()

9 MCMC sampling diagnostics

brms

bayesplot

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

See list of available diagnostics by name

available_mcmc()

bayesplot MCMC module:
  mcmc_acf
  mcmc_acf_bar
  mcmc_areas
  mcmc_areas_ridges
  mcmc_combo
  mcmc_dens
  mcmc_dens_chains
  mcmc_dens_overlay
  mcmc_hex
  mcmc_hist
  mcmc_hist_by_chain
  mcmc_intervals
  mcmc_neff
  mcmc_neff_hist
  mcmc_nuts_acceptance
  mcmc_nuts_divergence
  mcmc_nuts_energy
  mcmc_nuts_stepsize
  mcmc_nuts_treedepth
  mcmc_pairs
  mcmc_parcoord
  mcmc_rank_ecdf
  mcmc_rank_hist
  mcmc_rank_overlay
  mcmc_recover_hist
  mcmc_recover_intervals
  mcmc_recover_scatter
  mcmc_rhat
  mcmc_rhat_hist
  mcmc_scatter
  mcmc_trace
  mcmc_trace_highlight
  mcmc_violin

Of these, we will focus on:

• trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

pars <- hughes.brm3 %>% get_variables()
pars <- hughes.brm3 |>
  get_variables() |>
  str_subset("^b_.*|[sS]igma|^sd.*")
pars

 [1] "b_Intercept[1]"           "b_Intercept[2]"          
 [3] "b_Intercept[3]"           "b_Intercept[4]"          
 [5] "b_HABITATF"               "b_HABITATL"              
 [7] "b_HABITATU"               "b_SECTORCentral"         
 [9] "b_SECTORSouth"            "b_HABITATF:SECTORCentral"
[11] "b_HABITATL:SECTORCentral" "b_HABITATU:SECTORCentral"
[13] "b_HABITATF:SECTORSouth"   "b_HABITATL:SECTORSouth"  
[15] "b_HABITATU:SECTORSouth"   "sd_REEF__Intercept"      

hughes.brm3 |> mcmc_plot(type = "trace", variable = pars)

No divergences to plot.

The chains appear well mixed and very similar

• acf_bar (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

hughes.brm3 |> mcmc_plot(type = "acf_bar", variable = pars)

There is no evidence of auto-correlation in the MCMC samples

• rhat_hist: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

hughes.brm3 |> mcmc_plot(type = "rhat_hist")

Warning: Dropped 1 NAs from 'new_rhat(rhat)'.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

• neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

hughes.brm2 |> mcmc_plot(type = "neff_hist")

Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ratios all very high.

More diagnostics

hughes.brm3 |> mcmc_plot(type = "combo", variable = pars)

hughes.brm3 |> mcmc_plot(type = "violin", variable = pars)

stan plots

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

hughes.brm3 |> get_variables()

  [1] "b_Intercept[1]"            "b_Intercept[2]"           
  [3] "b_Intercept[3]"            "b_Intercept[4]"           
  [5] "b_HABITATF"                "b_HABITATL"               
  [7] "b_HABITATU"                "b_SECTORCentral"          
  [9] "b_SECTORSouth"             "b_HABITATF:SECTORCentral" 
 [11] "b_HABITATL:SECTORCentral"  "b_HABITATU:SECTORCentral" 
 [13] "b_HABITATF:SECTORSouth"    "b_HABITATL:SECTORSouth"   
 [15] "b_HABITATU:SECTORSouth"    "sd_REEF__Intercept"       
 [17] "disc"                      "Intercept[1]"             
 [19] "Intercept[2]"              "Intercept[3]"             
 [21] "Intercept[4]"              "r_REEF[09-304,Intercept]" 
 [23] "r_REEF[09-306,Intercept]"  "r_REEF[09-355,Intercept]" 
 [25] "r_REEF[09-357,Intercept]"  "r_REEF[09-358a,Intercept]"
 [27] "r_REEF[09-366,Intercept]"  "r_REEF[09-368,Intercept]" 
 [29] "r_REEF[09-372,Intercept]"  "r_REEF[09-376,Intercept]" 
 [31] "r_REEF[09-378,Intercept]"  "r_REEF[09-380,Intercept]" 
 [33] "r_REEF[09-394,Intercept]"  "r_REEF[09-406,Intercept]" 
 [35] "r_REEF[09-407,Intercept]"  "r_REEF[09-408,Intercept]" 
 [37] "r_REEF[09-409,Intercept]"  "r_REEF[09-410,Intercept]" 
 [39] "r_REEF[09-411,Intercept]"  "r_REEF[09-419a,Intercept]"
 [41] "r_REEF[09-423,Intercept]"  "r_REEF[09-424,Intercept]" 
 [43] "r_REEF[10-283b,Intercept]" "r_REEF[10-299,Intercept]" 
 [45] "r_REEF[10-302,Intercept]"  "r_REEF[10-303,Intercept]" 
 [47] "r_REEF[10-303d,Intercept]" "r_REEF[10-305,Intercept]" 
 [49] "r_REEF[10-318a,Intercept]" "r_REEF[10-318c,Intercept]"
 [51] "r_REEF[10-320,Intercept]"  "r_REEF[10-323,Intercept]" 
 [53] "r_REEF[10-325,Intercept]"  "r_REEF[10-326,Intercept]" 
 [55] "r_REEF[10-327,Intercept]"  "r_REEF[10-335,Intercept]" 
 [57] "r_REEF[10-336,Intercept]"  "r_REEF[10-337,Intercept]" 
 [59] "r_REEF[10-338,Intercept]"  "r_REEF[10-341,Intercept]" 
 [61] "r_REEF[10-347a,Intercept]" "r_REEF[10-349,Intercept]" 
 [63] "r_REEF[10-350,Intercept]"  "r_REEF[10-351,Intercept]" 
 [65] "r_REEF[10-352,Intercept]"  "r_REEF[10-353,Intercept]" 
 [67] "r_REEF[10-355,Intercept]"  "r_REEF[10-377,Intercept]" 
 [69] "r_REEF[10-381,Intercept]"  "r_REEF[10-384,Intercept]" 
 [71] "r_REEF[10-386,Intercept]"  "r_REEF[10-387,Intercept]" 
 [73] "r_REEF[10-390,Intercept]"  "r_REEF[10-391,Intercept]" 
 [75] "r_REEF[10-399,Intercept]"  "r_REEF[10-402,Intercept]" 
 [77] "r_REEF[10-406,Intercept]"  "r_REEF[10-412,Intercept]" 
 [79] "r_REEF[10-414,Intercept]"  "r_REEF[10-416b,Intercept]"
 [81] "r_REEF[10-418a,Intercept]" "r_REEF[10-419a,Intercept]"
 [83] "r_REEF[10-420,Intercept]"  "r_REEF[10-426,Intercept]" 
 [85] "r_REEF[10-441,Intercept]"  "r_REEF[10-449,Intercept]" 
 [87] "r_REEF[10-450,Intercept]"  "r_REEF[10-456a,Intercept]"
 [89] "r_REEF[10-458,Intercept]"  "r_REEF[10-458b,Intercept]"
 [91] "r_REEF[10-503,Intercept]"  "r_REEF[10-514,Intercept]" 
 [93] "r_REEF[10-518c,Intercept]" "r_REEF[10-521a,Intercept]"
 [95] "r_REEF[10-541,Intercept]"  "r_REEF[10-547b,Intercept]"
 [97] "r_REEF[10-562,Intercept]"  "r_REEF[10-565,Intercept]" 
 [99] "r_REEF[10-566a,Intercept]" "r_REEF[10-567,Intercept]" 
[101] "r_REEF[10-572,Intercept]"  "r_REEF[10-575,Intercept]" 
[103] "r_REEF[10-576,Intercept]"  "r_REEF[10-582a,Intercept]"
[105] "r_REEF[10-583,Intercept]"  "r_REEF[10-585,Intercept]" 
[107] "r_REEF[10-598a,Intercept]" "r_REEF[10-619,Intercept]" 
[109] "r_REEF[10-624,Intercept]"  "r_REEF[10-800,Intercept]" 
[111] "r_REEF[11-001,Intercept]"  "r_REEF[11-025,Intercept]" 
[113] "r_REEF[11-026,Intercept]"  "r_REEF[11-035,Intercept]" 
[115] "r_REEF[11-038,Intercept]"  "r_REEF[11-052,Intercept]" 
[117] "r_REEF[11-053,Intercept]"  "r_REEF[11-055,Intercept]" 
[119] "r_REEF[11-061,Intercept]"  "r_REEF[11-062,Intercept]" 
[121] "r_REEF[11-072,Intercept]"  "r_REEF[11-073,Intercept]" 
[123] "r_REEF[11-074,Intercept]"  "r_REEF[11-075,Intercept]" 
[125] "r_REEF[11-076,Intercept]"  "r_REEF[11-077,Intercept]" 
[127] "r_REEF[11-081,Intercept]"  "r_REEF[11-082,Intercept]" 
[129] "r_REEF[11-087b,Intercept]" "r_REEF[11-089b,Intercept]"
[131] "r_REEF[11-091,Intercept]"  "r_REEF[11-094,Intercept]" 
[133] "r_REEF[11-102,Intercept]"  "r_REEF[11-109a,Intercept]"
[135] "r_REEF[11-109b,Intercept]" "r_REEF[11-110,Intercept]" 
[137] "r_REEF[11-113,Intercept]"  "r_REEF[11-114,Intercept]" 
[139] "r_REEF[11-115,Intercept]"  "r_REEF[11-116,Intercept]" 
[141] "r_REEF[11-117,Intercept]"  "r_REEF[11-118a,Intercept]"
[143] "r_REEF[11-119,Intercept]"  "r_REEF[11-120,Intercept]" 
[145] "r_REEF[11-123,Intercept]"  "r_REEF[11-125,Intercept]" 
[147] "r_REEF[11-127b,Intercept]" "r_REEF[11-128,Intercept]" 
[149] "r_REEF[11-129,Intercept]"  "r_REEF[11-134,Intercept]" 
[151] "r_REEF[11-136,Intercept]"  "r_REEF[11-137,Intercept]" 
[153] "r_REEF[11-138,Intercept]"  "r_REEF[11-139,Intercept]" 
[155] "r_REEF[11-140,Intercept]"  "r_REEF[11-141,Intercept]" 
[157] "r_REEF[11-142,Intercept]"  "r_REEF[11-143,Intercept]" 
[159] "r_REEF[11-147,Intercept]"  "r_REEF[11-150,Intercept]" 
[161] "r_REEF[11-151,Intercept]"  "r_REEF[11-156,Intercept]" 
[163] "r_REEF[11-161,Intercept]"  "r_REEF[11-173,Intercept]" 
[165] "r_REEF[11-181,Intercept]"  "r_REEF[11-184a,Intercept]"
[167] "r_REEF[11-184b,Intercept]" "r_REEF[11-184c,Intercept]"
[169] "r_REEF[11-187,Intercept]"  "r_REEF[11-188,Intercept]" 
[171] "r_REEF[11-189,Intercept]"  "r_REEF[11-190a,Intercept]"
[173] "r_REEF[11-190b,Intercept]" "r_REEF[11-191,Intercept]" 
[175] "r_REEF[11-194,Intercept]"  "r_REEF[11-197,Intercept]" 
[177] "r_REEF[11-201,Intercept]"  "r_REEF[11-202,Intercept]" 
[179] "r_REEF[11-204,Intercept]"  "r_REEF[11-208a,Intercept]"
[181] "r_REEF[11-210a,Intercept]" "r_REEF[11-212,Intercept]" 
[183] "r_REEF[11-222a,Intercept]" "r_REEF[11-223,Intercept]" 
[185] "r_REEF[11-225,Intercept]"  "r_REEF[11-228,Intercept]" 
[187] "r_REEF[11-229a,Intercept]" "r_REEF[11-229b,Intercept]"
[189] "r_REEF[11-231,Intercept]"  "r_REEF[11-232,Intercept]" 
[191] "r_REEF[11-236,Intercept]"  "r_REEF[11-237a,Intercept]"
[193] "r_REEF[11-238,Intercept]"  "r_REEF[11-239,Intercept]" 
[195] "r_REEF[11-241,Intercept]"  "r_REEF[11-244c,Intercept]"
[197] "r_REEF[12-001,Intercept]"  "r_REEF[12-003a,Intercept]"
[199] "r_REEF[12-004,Intercept]"  "r_REEF[12-006,Intercept]" 
[201] "r_REEF[12-007,Intercept]"  "r_REEF[12-010,Intercept]" 
[203] "r_REEF[12-011,Intercept]"  "r_REEF[12-012,Intercept]" 
[205] "r_REEF[12-013,Intercept]"  "r_REEF[12-014,Intercept]" 
[207] "r_REEF[12-016a,Intercept]" "r_REEF[12-016b,Intercept]"
[209] "r_REEF[12-018,Intercept]"  "r_REEF[12-022,Intercept]" 
[211] "r_REEF[12-027a,Intercept]" "r_REEF[12-031,Intercept]" 
[213] "r_REEF[12-032,Intercept]"  "r_REEF[12-033,Intercept]" 
[215] "r_REEF[12-034,Intercept]"  "r_REEF[12-035,Intercept]" 
[217] "r_REEF[12-036,Intercept]"  "r_REEF[12-037,Intercept]" 
[219] "r_REEF[12-038a,Intercept]" "r_REEF[12-038c,Intercept]"
[221] "r_REEF[12-039,Intercept]"  "r_REEF[12-040,Intercept]" 
[223] "r_REEF[12-042,Intercept]"  "r_REEF[12-043a,Intercept]"
[225] "r_REEF[12-043b,Intercept]" "r_REEF[12-046,Intercept]" 
[227] "r_REEF[12-048,Intercept]"  "r_REEF[12-051,Intercept]" 
[229] "r_REEF[12-053,Intercept]"  "r_REEF[12-056,Intercept]" 
[231] "r_REEF[12-058,Intercept]"  "r_REEF[12-059,Intercept]" 
[233] "r_REEF[12-061,Intercept]"  "r_REEF[12-068,Intercept]" 
[235] "r_REEF[12-069,Intercept]"  "r_REEF[12-070,Intercept]" 
[237] "r_REEF[12-071,Intercept]"  "r_REEF[12-075,Intercept]" 
[239] "r_REEF[12-076,Intercept]"  "r_REEF[12-077,Intercept]" 
[241] "r_REEF[12-078,Intercept]"  "r_REEF[12-094,Intercept]" 
[243] "r_REEF[12-096,Intercept]"  "r_REEF[12-098,Intercept]" 
[245] "r_REEF[12-100,Intercept]"  "r_REEF[12-101,Intercept]" 
[247] "r_REEF[12-102,Intercept]"  "r_REEF[12-104,Intercept]" 
[249] "r_REEF[12-105,Intercept]"  "r_REEF[12-107,Intercept]" 
[251] "r_REEF[12-110,Intercept]"  "r_REEF[12-111,Intercept]" 
[253] "r_REEF[12-112,Intercept]"  "r_REEF[12-114,Intercept]" 
[255] "r_REEF[12-115,Intercept]"  "r_REEF[12-116,Intercept]" 
[257] "r_REEF[12-118,Intercept]"  "r_REEF[12-119,Intercept]" 
[259] "r_REEF[12-127,Intercept]"  "r_REEF[12-129,Intercept]" 
[261] "r_REEF[12-131,Intercept]"  "r_REEF[12-134,Intercept]" 
[263] "r_REEF[12-137,Intercept]"  "r_REEF[12-143,Intercept]" 
[265] "r_REEF[12-145,Intercept]"  "r_REEF[13-002,Intercept]" 
[267] "r_REEF[13-005,Intercept]"  "r_REEF[13-006,Intercept]" 
[269] "r_REEF[13-007,Intercept]"  "r_REEF[13-012,Intercept]" 
[271] "r_REEF[13-014,Intercept]"  "r_REEF[13-019,Intercept]" 
[273] "r_REEF[13-022,Intercept]"  "r_REEF[13-028,Intercept]" 
[275] "r_REEF[13-029,Intercept]"  "r_REEF[13-031,Intercept]" 
[277] "r_REEF[13-032,Intercept]"  "r_REEF[13-034,Intercept]" 
[279] "r_REEF[13-040,Intercept]"  "r_REEF[13-041,Intercept]" 
[281] "r_REEF[13-045,Intercept]"  "r_REEF[13-050a,Intercept]"
[283] "r_REEF[13-055,Intercept]"  "r_REEF[13-056,Intercept]" 
[285] "r_REEF[13-060,Intercept]"  "r_REEF[13-061a,Intercept]"
[287] "r_REEF[13-061b,Intercept]" "r_REEF[13-063,Intercept]" 
[289] "r_REEF[13-069,Intercept]"  "r_REEF[13-072,Intercept]" 
[291] "r_REEF[13-074,Intercept]"  "r_REEF[13-075,Intercept]" 
[293] "r_REEF[13-076,Intercept]"  "r_REEF[13-077,Intercept]" 
[295] "r_REEF[13-078,Intercept]"  "r_REEF[13-079c,Intercept]"
[297] "r_REEF[13-081,Intercept]"  "r_REEF[13-083,Intercept]" 
[299] "r_REEF[13-087,Intercept]"  "r_REEF[13-088,Intercept]" 
[301] "r_REEF[13-091,Intercept]"  "r_REEF[13-093b,Intercept]"
[303] "r_REEF[13-093c,Intercept]" "r_REEF[13-097,Intercept]" 
[305] "r_REEF[13-108,Intercept]"  "r_REEF[13-111,Intercept]" 
[307] "r_REEF[13-116,Intercept]"  "r_REEF[13-117b,Intercept]"
[309] "r_REEF[13-118a,Intercept]" "r_REEF[13-119,Intercept]" 
[311] "r_REEF[13-120,Intercept]"  "r_REEF[13-121,Intercept]" 
[313] "r_REEF[13-122,Intercept]"  "r_REEF[13-125,Intercept]" 
[315] "r_REEF[13-127,Intercept]"  "r_REEF[13-129,Intercept]" 
[317] "r_REEF[13-130,Intercept]"  "r_REEF[13-131,Intercept]" 
[319] "r_REEF[13-133,Intercept]"  "r_REEF[14-003,Intercept]" 
[321] "r_REEF[14-006,Intercept]"  "r_REEF[14-008,Intercept]" 
[323] "r_REEF[14-012,Intercept]"  "r_REEF[14-013,Intercept]" 
[325] "r_REEF[14-016,Intercept]"  "r_REEF[14-017,Intercept]" 
[327] "r_REEF[14-022,Intercept]"  "r_REEF[14-026a,Intercept]"
[329] "r_REEF[14-026b,Intercept]" "r_REEF[14-028,Intercept]" 
[331] "r_REEF[14-029,Intercept]"  "r_REEF[14-032,Intercept]" 
[333] "r_REEF[14-033a,Intercept]" "r_REEF[14-034,Intercept]" 
[335] "r_REEF[14-038,Intercept]"  "r_REEF[14-039,Intercept]" 
[337] "r_REEF[14-042a,Intercept]" "r_REEF[14-045,Intercept]" 
[339] "r_REEF[14-047,Intercept]"  "r_REEF[14-049,Intercept]" 
[341] "r_REEF[14-051,Intercept]"  "r_REEF[14-052,Intercept]" 
[343] "r_REEF[14-053,Intercept]"  "r_REEF[14-055,Intercept]" 
[345] "r_REEF[14-061,Intercept]"  "r_REEF[14-063,Intercept]" 
[347] "r_REEF[14-064,Intercept]"  "r_REEF[14-065,Intercept]" 
[349] "r_REEF[14-066,Intercept]"  "r_REEF[14-068,Intercept]" 
[351] "r_REEF[14-070,Intercept]"  "r_REEF[14-071a,Intercept]"
[353] "r_REEF[14-072,Intercept]"  "r_REEF[14-073,Intercept]" 
[355] "r_REEF[14-074,Intercept]"  "r_REEF[14-075,Intercept]" 
[357] "r_REEF[14-077a,Intercept]" "r_REEF[14-078,Intercept]" 
[359] "r_REEF[14-079,Intercept]"  "r_REEF[14-082,Intercept]" 
[361] "r_REEF[14-085,Intercept]"  "r_REEF[14-086,Intercept]" 
[363] "r_REEF[14-088,Intercept]"  "r_REEF[14-089,Intercept]" 
[365] "r_REEF[14-090,Intercept]"  "r_REEF[14-091,Intercept]" 
[367] "r_REEF[14-093,Intercept]"  "r_REEF[14-094,Intercept]" 
[369] "r_REEF[14-096,Intercept]"  "r_REEF[14-097,Intercept]" 
[371] "r_REEF[14-099,Intercept]"  "r_REEF[14-101,Intercept]" 
[373] "r_REEF[14-102,Intercept]"  "r_REEF[14-103,Intercept]" 
[375] "r_REEF[14-104,Intercept]"  "r_REEF[14-106,Intercept]" 
[377] "r_REEF[14-116a,Intercept]" "r_REEF[14-116b,Intercept]"
[379] "r_REEF[14-116c,Intercept]" "r_REEF[14-116d,Intercept]"
[381] "r_REEF[14-120d,Intercept]" "r_REEF[14-122a,Intercept]"
[383] "r_REEF[14-123,Intercept]"  "r_REEF[14-126,Intercept]" 
[385] "r_REEF[14-132a,Intercept]" "r_REEF[14-132b,Intercept]"
[387] "r_REEF[14-133,Intercept]"  "r_REEF[14-135,Intercept]" 
[389] "r_REEF[14-137,Intercept]"  "r_REEF[14-138,Intercept]" 
[391] "r_REEF[14-139,Intercept]"  "r_REEF[14-140,Intercept]" 
[393] "r_REEF[14-143,Intercept]"  "r_REEF[14-146,Intercept]" 
[395] "r_REEF[14-147a,Intercept]" "r_REEF[14-147b,Intercept]"
[397] "r_REEF[14-150,Intercept]"  "r_REEF[14-154,Intercept]" 
[399] "r_REEF[15-002,Intercept]"  "r_REEF[15-005,Intercept]" 
[401] "r_REEF[15-006,Intercept]"  "r_REEF[15-008,Intercept]" 
[403] "r_REEF[15-009,Intercept]"  "r_REEF[15-012,Intercept]" 
[405] "r_REEF[15-013,Intercept]"  "r_REEF[15-017a,Intercept]"
[407] "r_REEF[15-018,Intercept]"  "r_REEF[15-019,Intercept]" 
[409] "r_REEF[15-021a,Intercept]" "r_REEF[15-022,Intercept]" 
[411] "r_REEF[15-024,Intercept]"  "r_REEF[15-026,Intercept]" 
[413] "r_REEF[15-028,Intercept]"  "r_REEF[15-032,Intercept]" 
[415] "r_REEF[15-033,Intercept]"  "r_REEF[15-034,Intercept]" 
[417] "r_REEF[15-038,Intercept]"  "r_REEF[15-039,Intercept]" 
[419] "r_REEF[15-043,Intercept]"  "r_REEF[15-046,Intercept]" 
[421] "r_REEF[15-049,Intercept]"  "r_REEF[15-050,Intercept]" 
[423] "r_REEF[15-052,Intercept]"  "r_REEF[15-063,Intercept]" 
[425] "r_REEF[15-064,Intercept]"  "r_REEF[15-065,Intercept]" 
[427] "r_REEF[15-066,Intercept]"  "r_REEF[15-071a,Intercept]"
[429] "r_REEF[15-075a,Intercept]" "r_REEF[15-078,Intercept]" 
[431] "r_REEF[15-080,Intercept]"  "r_REEF[15-085,Intercept]" 
[433] "r_REEF[15-086,Intercept]"  "r_REEF[15-087,Intercept]" 
[435] "r_REEF[15-088,Intercept]"  "r_REEF[15-089,Intercept]" 
[437] "r_REEF[15-090,Intercept]"  "r_REEF[15-092,Intercept]" 
[439] "r_REEF[15-093,Intercept]"  "r_REEF[15-094,Intercept]" 
[441] "r_REEF[15-095,Intercept]"  "r_REEF[15-096,Intercept]" 
[443] "r_REEF[15-098,Intercept]"  "r_REEF[15-099a,Intercept]"
[445] "r_REEF[15-099d,Intercept]" "r_REEF[16-015,Intercept]" 
[447] "r_REEF[16-019,Intercept]"  "r_REEF[16-020a,Intercept]"
[449] "r_REEF[16-023,Intercept]"  "r_REEF[16-025,Intercept]" 
[451] "r_REEF[16-026,Intercept]"  "r_REEF[16-029,Intercept]" 
[453] "r_REEF[16-030,Intercept]"  "r_REEF[16-032,Intercept]" 
[455] "r_REEF[16-040,Intercept]"  "r_REEF[16-043a,Intercept]"
[457] "r_REEF[16-046,Intercept]"  "r_REEF[16-049,Intercept]" 
[459] "r_REEF[16-054b,Intercept]" "r_REEF[16-054c,Intercept]"
[461] "r_REEF[16-057,Intercept]"  "r_REEF[16-060,Intercept]" 
[463] "r_REEF[16-068,Intercept]"  "r_REEF[16-071,Intercept]" 
[465] "r_REEF[16-073,Intercept]"  "r_REEF[17-001a,Intercept]"
[467] "r_REEF[17-004,Intercept]"  "r_REEF[17-005,Intercept]" 
[469] "r_REEF[17-006,Intercept]"  "r_REEF[17-008,Intercept]" 
[471] "r_REEF[17-010,Intercept]"  "r_REEF[17-011,Intercept]" 
[473] "r_REEF[17-014,Intercept]"  "r_REEF[17-016a,Intercept]"
[475] "r_REEF[17-017,Intercept]"  "r_REEF[17-018a,Intercept]"
[477] "r_REEF[17-023,Intercept]"  "r_REEF[17-024,Intercept]" 
[479] "r_REEF[17-034,Intercept]"  "r_REEF[17-035,Intercept]" 
[481] "r_REEF[17-037,Intercept]"  "r_REEF[17-042,Intercept]" 
[483] "r_REEF[17-044,Intercept]"  "r_REEF[17-047,Intercept]" 
[485] "r_REEF[17-051,Intercept]"  "r_REEF[17-057,Intercept]" 
[487] "r_REEF[17-059a,Intercept]" "r_REEF[17-062,Intercept]" 
[489] "r_REEF[17-063a,Intercept]" "r_REEF[17-064,Intercept]" 
[491] "r_REEF[17-065,Intercept]"  "r_REEF[17-066,Intercept]" 
[493] "r_REEF[17-067,Intercept]"  "r_REEF[17-068,Intercept]" 
[495] "r_REEF[17-069,Intercept]"  "r_REEF[18-017,Intercept]" 
[497] "r_REEF[18-018,Intercept]"  "r_REEF[18-019,Intercept]" 
[499] "r_REEF[18-023,Intercept]"  "r_REEF[18-024,Intercept]" 
[501] "r_REEF[18-026,Intercept]"  "r_REEF[18-027,Intercept]" 
[503] "r_REEF[18-029,Intercept]"  "r_REEF[18-030,Intercept]" 
[505] "r_REEF[18-031,Intercept]"  "r_REEF[18-032,Intercept]" 
[507] "r_REEF[18-033,Intercept]"  "r_REEF[18-034,Intercept]" 
[509] "r_REEF[18-037,Intercept]"  "r_REEF[18-039,Intercept]" 
[511] "r_REEF[18-042,Intercept]"  "r_REEF[18-043,Intercept]" 
[513] "r_REEF[18-048,Intercept]"  "r_REEF[18-049a,Intercept]"
[515] "r_REEF[18-049b,Intercept]" "r_REEF[18-049c,Intercept]"
[517] "r_REEF[18-049e,Intercept]" "r_REEF[18-054b,Intercept]"
[519] "r_REEF[18-054d,Intercept]" "r_REEF[18-054e,Intercept]"
[521] "r_REEF[18-054f,Intercept]" "r_REEF[18-065,Intercept]" 
[523] "r_REEF[18-069a,Intercept]" "r_REEF[18-069b,Intercept]"
[525] "r_REEF[18-070,Intercept]"  "r_REEF[18-071,Intercept]" 
[527] "r_REEF[18-075,Intercept]"  "r_REEF[18-077,Intercept]" 
[529] "r_REEF[18-078,Intercept]"  "r_REEF[18-084a,Intercept]"
[531] "r_REEF[18-084b,Intercept]" "r_REEF[18-085a,Intercept]"
[533] "r_REEF[18-086,Intercept]"  "r_REEF[18-091,Intercept]" 
[535] "r_REEF[18-096,Intercept]"  "r_REEF[18-099,Intercept]" 
[537] "r_REEF[18-100a,Intercept]" "r_REEF[18-100b,Intercept]"
[539] "r_REEF[18-101,Intercept]"  "r_REEF[18-103,Intercept]" 
[541] "r_REEF[18-105a,Intercept]" "r_REEF[18-105b,Intercept]"
[543] "r_REEF[18-106,Intercept]"  "r_REEF[18-115,Intercept]" 
[545] "r_REEF[18-118,Intercept]"  "r_REEF[18-120,Intercept]" 
[547] "r_REEF[19-009a,Intercept]" "r_REEF[19-009g,Intercept]"
[549] "r_REEF[19-009j,Intercept]" "r_REEF[19-019,Intercept]" 
[551] "r_REEF[19-024,Intercept]"  "r_REEF[19-038a,Intercept]"
[553] "r_REEF[19-038b,Intercept]" "r_REEF[19-038f,Intercept]"
[555] "r_REEF[19-043,Intercept]"  "r_REEF[19-045,Intercept]" 
[557] "r_REEF[19-048,Intercept]"  "r_REEF[19-049,Intercept]" 
[559] "r_REEF[19-050,Intercept]"  "r_REEF[19-051,Intercept]" 
[561] "r_REEF[19-054,Intercept]"  "r_REEF[19-059,Intercept]" 
[563] "r_REEF[19-061,Intercept]"  "r_REEF[19-063a,Intercept]"
[565] "r_REEF[19-065,Intercept]"  "r_REEF[19-066,Intercept]" 
[567] "r_REEF[19-067,Intercept]"  "r_REEF[19-071,Intercept]" 
[569] "r_REEF[19-072a,Intercept]" "r_REEF[19-072b,Intercept]"
[571] "r_REEF[19-072e,Intercept]" "r_REEF[19-074b,Intercept]"
[573] "r_REEF[19-076,Intercept]"  "r_REEF[19-081,Intercept]" 
[575] "r_REEF[19-082,Intercept]"  "r_REEF[19-091,Intercept]" 
[577] "r_REEF[19-096,Intercept]"  "r_REEF[19-097,Intercept]" 
[579] "r_REEF[19-098,Intercept]"  "r_REEF[19-099,Intercept]" 
[581] "r_REEF[19-107,Intercept]"  "r_REEF[19-109a,Intercept]"
[583] "r_REEF[19-113,Intercept]"  "r_REEF[19-114,Intercept]" 
[585] "r_REEF[19-116,Intercept]"  "r_REEF[19-123,Intercept]" 
[587] "r_REEF[19-127,Intercept]"  "r_REEF[19-128,Intercept]" 
[589] "r_REEF[19-129,Intercept]"  "r_REEF[19-134,Intercept]" 
[591] "r_REEF[19-135,Intercept]"  "r_REEF[19-136a,Intercept]"
[593] "r_REEF[19-139,Intercept]"  "r_REEF[19-141,Intercept]" 
[595] "r_REEF[19-146,Intercept]"  "r_REEF[19-147,Intercept]" 
[597] "r_REEF[19-151,Intercept]"  "r_REEF[19-152,Intercept]" 
[599] "r_REEF[19-157,Intercept]"  "r_REEF[19-159,Intercept]" 
[601] "r_REEF[19-167,Intercept]"  "r_REEF[19-168,Intercept]" 
[603] "r_REEF[19-169,Intercept]"  "r_REEF[19-187,Intercept]" 
[605] "r_REEF[19-197,Intercept]"  "r_REEF[19-204,Intercept]" 
[607] "r_REEF[19-206,Intercept]"  "r_REEF[19-207,Intercept]" 
[609] "r_REEF[19-208,Intercept]"  "r_REEF[19-209,Intercept]" 
[611] "r_REEF[19-212,Intercept]"  "r_REEF[19-214,Intercept]" 
[613] "r_REEF[19-216,Intercept]"  "r_REEF[19-220,Intercept]" 
[615] "r_REEF[19-221,Intercept]"  "r_REEF[19-222,Intercept]" 
[617] "r_REEF[19-225,Intercept]"  "r_REEF[19-227,Intercept]" 
[619] "r_REEF[19-246,Intercept]"  "r_REEF[20-014,Intercept]" 
[621] "r_REEF[20-017,Intercept]"  "r_REEF[20-028h,Intercept]"
[623] "r_REEF[20-078a,Intercept]" "r_REEF[20-104,Intercept]" 
[625] "r_REEF[20-105,Intercept]"  "r_REEF[20-112,Intercept]" 
[627] "r_REEF[20-115,Intercept]"  "r_REEF[20-119,Intercept]" 
[629] "r_REEF[20-121,Intercept]"  "r_REEF[20-125,Intercept]" 
[631] "r_REEF[20-126,Intercept]"  "r_REEF[20-127,Intercept]" 
[633] "r_REEF[20-128,Intercept]"  "r_REEF[20-129,Intercept]" 
[635] "r_REEF[20-131,Intercept]"  "r_REEF[20-132,Intercept]" 
[637] "r_REEF[20-134,Intercept]"  "r_REEF[20-135,Intercept]" 
[639] "r_REEF[20-136,Intercept]"  "r_REEF[20-137,Intercept]" 
[641] "r_REEF[20-138,Intercept]"  "r_REEF[20-167,Intercept]" 
[643] "r_REEF[20-172,Intercept]"  "r_REEF[20-179,Intercept]" 
[645] "r_REEF[20-182,Intercept]"  "r_REEF[20-185a,Intercept]"
[647] "r_REEF[20-185c,Intercept]" "r_REEF[20-186,Intercept]" 
[649] "r_REEF[20-187,Intercept]"  "r_REEF[20-194a,Intercept]"
[651] "r_REEF[20-197,Intercept]"  "r_REEF[20-269,Intercept]" 
[653] "r_REEF[20-291,Intercept]"  "r_REEF[20-292,Intercept]" 
[655] "r_REEF[20-295,Intercept]"  "r_REEF[20-296,Intercept]" 
[657] "r_REEF[20-298,Intercept]"  "r_REEF[20-299a,Intercept]"
[659] "r_REEF[20-299c,Intercept]" "r_REEF[20-301a,Intercept]"
[661] "r_REEF[20-301b,Intercept]" "r_REEF[20-302,Intercept]" 
[663] "r_REEF[20-304,Intercept]"  "r_REEF[20-305,Intercept]" 
[665] "r_REEF[20-306,Intercept]"  "r_REEF[20-307,Intercept]" 
[667] "r_REEF[20-308,Intercept]"  "r_REEF[20-322,Intercept]" 
[669] "r_REEF[20-324,Intercept]"  "r_REEF[20-326,Intercept]" 
[671] "r_REEF[20-335,Intercept]"  "r_REEF[20-339a,Intercept]"
[673] "r_REEF[20-339b,Intercept]" "r_REEF[20-348,Intercept]" 
[675] "r_REEF[20-353,Intercept]"  "r_REEF[20-369a,Intercept]"
[677] "r_REEF[20-370a,Intercept]" "r_REEF[20-370b,Intercept]"
[679] "r_REEF[20-377,Intercept]"  "r_REEF[20-378a,Intercept]"
[681] "r_REEF[20-380,Intercept]"  "r_REEF[20-382,Intercept]" 
[683] "r_REEF[20-383,Intercept]"  "r_REEF[20-385,Intercept]" 
[685] "r_REEF[20-386,Intercept]"  "r_REEF[20-388a,Intercept]"
[687] "r_REEF[20-389,Intercept]"  "r_REEF[20-391,Intercept]" 
[689] "r_REEF[20-469,Intercept]"  "r_REEF[20-470,Intercept]" 
[691] "r_REEF[20-474,Intercept]"  "r_REEF[20-568,Intercept]" 
[693] "r_REEF[21-025,Intercept]"  "r_REEF[21-050,Intercept]" 
[695] "r_REEF[21-054,Intercept]"  "r_REEF[21-056b,Intercept]"
[697] "r_REEF[21-057,Intercept]"  "r_REEF[21-061,Intercept]" 
[699] "r_REEF[21-062,Intercept]"  "r_REEF[21-064,Intercept]" 
[701] "r_REEF[21-065,Intercept]"  "r_REEF[21-068,Intercept]" 
[703] "r_REEF[21-070,Intercept]"  "r_REEF[21-071,Intercept]" 
[705] "r_REEF[21-072,Intercept]"  "r_REEF[21-074,Intercept]" 
[707] "r_REEF[21-076,Intercept]"  "r_REEF[21-081,Intercept]" 
[709] "r_REEF[21-082,Intercept]"  "r_REEF[21-086,Intercept]" 
[711] "r_REEF[21-087,Intercept]"  "r_REEF[21-091,Intercept]" 
[713] "r_REEF[21-092,Intercept]"  "r_REEF[21-093,Intercept]" 
[715] "r_REEF[21-094,Intercept]"  "r_REEF[21-100,Intercept]" 
[717] "r_REEF[21-101,Intercept]"  "r_REEF[21-109,Intercept]" 
[719] "r_REEF[21-115,Intercept]"  "r_REEF[21-116,Intercept]" 
[721] "r_REEF[21-128,Intercept]"  "r_REEF[21-129,Intercept]" 
[723] "r_REEF[21-131,Intercept]"  "r_REEF[21-132,Intercept]" 
[725] "r_REEF[21-133,Intercept]"  "r_REEF[21-134,Intercept]" 
[727] "r_REEF[21-135,Intercept]"  "r_REEF[21-136,Intercept]" 
[729] "r_REEF[21-139,Intercept]"  "r_REEF[21-143,Intercept]" 
[731] "r_REEF[21-144,Intercept]"  "r_REEF[21-145a,Intercept]"
[733] "r_REEF[21-146,Intercept]"  "r_REEF[21-147,Intercept]" 
[735] "r_REEF[21-148,Intercept]"  "r_REEF[21-149,Intercept]" 
[737] "r_REEF[21-150,Intercept]"  "r_REEF[21-151,Intercept]" 
[739] "r_REEF[21-152,Intercept]"  "r_REEF[21-157,Intercept]" 
[741] "r_REEF[21-158,Intercept]"  "r_REEF[21-159,Intercept]" 
[743] "r_REEF[21-161,Intercept]"  "r_REEF[21-165,Intercept]" 
[745] "r_REEF[21-172,Intercept]"  "r_REEF[21-180,Intercept]" 
[747] "r_REEF[21-181,Intercept]"  "r_REEF[21-185,Intercept]" 
[749] "r_REEF[21-236,Intercept]"  "r_REEF[21-237,Intercept]" 
[751] "r_REEF[21-238,Intercept]"  "r_REEF[21-241,Intercept]" 
[753] "r_REEF[21-242a,Intercept]" "r_REEF[21-247,Intercept]" 
[755] "r_REEF[21-258,Intercept]"  "r_REEF[21-263,Intercept]" 
[757] "r_REEF[21-289a,Intercept]" "r_REEF[21-290,Intercept]" 
[759] "r_REEF[21-294,Intercept]"  "r_REEF[21-295,Intercept]" 
[761] "r_REEF[21-296,Intercept]"  "r_REEF[21-299,Intercept]" 
[763] "r_REEF[21-300a,Intercept]" "r_REEF[21-303,Intercept]" 
[765] "r_REEF[21-304,Intercept]"  "r_REEF[21-306,Intercept]" 
[767] "r_REEF[21-334,Intercept]"  "r_REEF[21-386,Intercept]" 
[769] "r_REEF[21-394,Intercept]"  "r_REEF[21-403b,Intercept]"
[771] "r_REEF[21-433,Intercept]"  "r_REEF[21-449,Intercept]" 
[773] "r_REEF[21-450,Intercept]"  "r_REEF[21-453,Intercept]" 
[775] "r_REEF[21-460,Intercept]"  "r_REEF[21-464,Intercept]" 
[777] "r_REEF[21-465,Intercept]"  "r_REEF[21-466,Intercept]" 
[779] "r_REEF[21-467,Intercept]"  "r_REEF[21-470,Intercept]" 
[781] "r_REEF[21-474,Intercept]"  "r_REEF[21-479,Intercept]" 
[783] "r_REEF[21-491,Intercept]"  "r_REEF[21-496,Intercept]" 
[785] "r_REEF[21-526,Intercept]"  "r_REEF[21-529,Intercept]" 
[787] "r_REEF[21-530,Intercept]"  "r_REEF[21-531,Intercept]" 
[789] "r_REEF[21-532,Intercept]"  "r_REEF[21-538,Intercept]" 
[791] "r_REEF[21-540,Intercept]"  "r_REEF[21-544,Intercept]" 
[793] "r_REEF[21-546,Intercept]"  "r_REEF[21-551,Intercept]" 
[795] "r_REEF[21-556,Intercept]"  "r_REEF[21-560,Intercept]" 
[797] "r_REEF[21-562,Intercept]"  "r_REEF[21-565,Intercept]" 
[799] "r_REEF[21-566,Intercept]"  "r_REEF[21-570,Intercept]" 
[801] "r_REEF[21-575,Intercept]"  "r_REEF[21-579,Intercept]" 
[803] "r_REEF[21-582,Intercept]"  "r_REEF[21-586,Intercept]" 
[805] "r_REEF[21-591,Intercept]"  "r_REEF[21-761,Intercept]" 
[807] "r_REEF[22-084,Intercept]"  "r_REEF[22-086,Intercept]" 
[809] "r_REEF[22-090,Intercept]"  "r_REEF[22-092,Intercept]" 
[811] "r_REEF[22-093,Intercept]"  "r_REEF[22-103,Intercept]" 
[813] "r_REEF[22-106,Intercept]"  "r_REEF[22-115,Intercept]" 
[815] "r_REEF[22-116,Intercept]"  "r_REEF[22-121,Intercept]" 
[817] "r_REEF[22-122,Intercept]"  "r_REEF[22-123,Intercept]" 
[819] "r_REEF[22-124,Intercept]"  "r_REEF[22-125,Intercept]" 
[821] "r_REEF[22-129,Intercept]"  "r_REEF[22-131,Intercept]" 
[823] "r_REEF[22-133,Intercept]"  "r_REEF[22-134,Intercept]" 
[825] "r_REEF[22-139,Intercept]"  "r_REEF[22-140,Intercept]" 
[827] "r_REEF[22-141,Intercept]"  "r_REEF[22-142,Intercept]" 
[829] "r_REEF[22-143,Intercept]"  "r_REEF[22-144,Intercept]" 
[831] "prior_Intercept"           "prior_b"                  
[833] "prior_sd_REEF"             "lprior"                   
[835] "lp__"                      "accept_stat__"            
[837] "stepsize__"                "treedepth__"              
[839] "n_leapfrog__"              "divergent__"              
[841] "energy__"                 

pars <- hughes.brm3 |>
  get_variables() |>
  str_subset("^b_.*|[sS]igma|^sd.*")
pars

 [1] "b_Intercept[1]"           "b_Intercept[2]"          
 [3] "b_Intercept[3]"           "b_Intercept[4]"          
 [5] "b_HABITATF"               "b_HABITATL"              
 [7] "b_HABITATU"               "b_SECTORCentral"         
 [9] "b_SECTORSouth"            "b_HABITATF:SECTORCentral"
[11] "b_HABITATL:SECTORCentral" "b_HABITATU:SECTORCentral"
[13] "b_HABITATF:SECTORSouth"   "b_HABITATL:SECTORSouth"  
[15] "b_HABITATU:SECTORSouth"   "sd_REEF__Intercept"      

hughes.brm3$fit |>
  stan_trace(pars = pars)

The chains appear well mixed and very similar

• stan_acf (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain

hughes.brm3$fit |>
  stan_ac(pars = pars)

There is no evidence of auto-correlation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

hughes.brm3$fit |> stan_rhat()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_bin()`).

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

hughes.brm3$fit |> stan_ess()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 1 row containing non-finite outside the scale range
(`stat_bin()`).

Ratios all very high.

hughes.brm3$fit |>
  stan_dens(separate_chains = TRUE, pars = pars)

10 Model validation

brms

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.

See list of available diagnostics by name

available_ppc()

bayesplot PPC module:
  ppc_bars
  ppc_bars_grouped
  ppc_boxplot
  ppc_dens
  ppc_dens_overlay
  ppc_dens_overlay_grouped
  ppc_dots
  ppc_ecdf_overlay
  ppc_ecdf_overlay_grouped
  ppc_error_binned
  ppc_error_hist
  ppc_error_hist_grouped
  ppc_error_scatter
  ppc_error_scatter_avg
  ppc_error_scatter_avg_grouped
  ppc_error_scatter_avg_vs_x
  ppc_freqpoly
  ppc_freqpoly_grouped
  ppc_hist
  ppc_intervals
  ppc_intervals_grouped
  ppc_km_overlay
  ppc_km_overlay_grouped
  ppc_loo_intervals
  ppc_loo_pit_ecdf
  ppc_loo_pit_overlay
  ppc_loo_pit_qq
  ppc_loo_ribbon
  ppc_pit_ecdf
  ppc_pit_ecdf_grouped
  ppc_ribbon
  ppc_ribbon_grouped
  ppc_rootogram
  ppc_scatter
  ppc_scatter_avg
  ppc_scatter_avg_grouped
  ppc_stat
  ppc_stat_2d
  ppc_stat_freqpoly
  ppc_stat_freqpoly_grouped
  ppc_stat_grouped
  ppc_violin_grouped

• 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.

hughes.brm3 |> pp_check(type = "dens_overlay", ndraws = 100)

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

## hughes.brm3 %>% pp_check(type = 'error_scatter_avg')

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.

hughes.brm3 %>% pp_check(group = "REEF", type = "intervals")

Using all posterior draws for ppc type 'intervals' by default.

Warning: The following arguments were unrecognized and ignored: group

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

# library(shinystan)
# launch_shinystan(hughes.brm2)

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

hughes.resids <- make_brms_dharma_res(hughes.brm3, integerResponse = FALSE)
wrap_elements(~ testUniformity(hughes.resids)) +
  wrap_elements(~ plotResiduals(hughes.resids, form = factor(rep(1, nrow(hughes))))) +
  wrap_elements(~ plotResiduals(hughes.resids, quantreg = TRUE)) +
  wrap_elements(~ plotResiduals(hughes.resids, quantreg = TRUE, asFactor = FALSE)) +
  wrap_elements(~ testDispersion(hughes.resids))

11 Model investigation

brms

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

hughes.brm3 |>
  conditional_effects("HABITAT",
    conditions = make_conditions(hughes.brm3, "SECTOR"),
    categorical = TRUE
  )

hughes.brm3 %>% summary()

 Family: cumulative 
  Links: mu = logit; disc = identity 
Formula: oSCORE ~ HABITAT * SECTOR + (1 | REEF) 
   Data: hughes (Number of observations: 3394) 
  Draws: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
         total post-warmup draws = 2400

Multilevel Hyperparameters:
~REEF (Number of levels: 809) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.71      0.11     2.49     2.94 1.00     1963     2291

Regression Coefficients:
                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept[1]              -8.49      0.24    -8.96    -8.03 1.00     2104
Intercept[2]              -5.34      0.20    -5.75    -4.95 1.00     2161
Intercept[3]              -2.74      0.17    -3.07    -2.41 1.00     2074
Intercept[4]              -0.40      0.16    -0.71    -0.09 1.00     2118
HABITATF                  -0.34      0.18    -0.70     0.01 1.00     2262
HABITATL                  -1.34      0.18    -1.69    -1.00 1.00     2257
HABITATU                   1.18      0.35     0.49     1.90 1.00     2416
SECTORCentral             -1.72      0.33    -2.35    -1.10 1.00     1964
SECTORSouth               -6.46      0.26    -7.00    -5.94 1.00     1868
HABITATF:SECTORCentral    -0.64      0.31    -1.25    -0.03 1.00     2391
HABITATL:SECTORCentral     0.22      0.31    -0.39     0.84 1.00     2427
HABITATU:SECTORCentral    -0.59      0.51    -1.58     0.40 1.00     2455
HABITATF:SECTORSouth      -0.40      0.23    -0.86     0.06 1.00     2177
HABITATL:SECTORSouth      -0.50      0.23    -0.98    -0.05 1.00     2295
HABITATU:SECTORSouth      -0.36      0.40    -1.17     0.42 1.00     2428
                       Tail_ESS
Intercept[1]               2217
Intercept[2]               2226
Intercept[3]               2143
Intercept[4]               2307
HABITATF                   2289
HABITATL                   2157
HABITATU                   2246
SECTORCentral              1796
SECTORSouth                1990
HABITATF:SECTORCentral     2252
HABITATL:SECTORCentral     2250
HABITATU:SECTORCentral     2339
HABITATF:SECTORSouth       2259
HABITATL:SECTORSouth       2075
HABITATU:SECTORSouth       2454

Further Distributional Parameters:
     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
disc     1.00      0.00     1.00     1.00   NA       NA       NA

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

• -8.49 log odds of being in the first category vs a higher category (0.0002)
• -5.35 (0.0047)
• -2.75 (0.06)
• -0.40 (0.67) odds of being in the 4th category (3) vs higher
• -0.35 (0.7) odds of being in a higher category are 30% lower in the F habitat than the C habitat (In the North)
• -1.34 (0.26) odds of being in a higher category are 64% lower in the L habitat than the C habitat (in the North)
• 1.20 (3.3) odds of being in a higher category are 232% higher in the U habitat than the C habitat (in the North)
• -1.72 (0.179) odds of being in a higher category are 82% lower in central vs northern (for C habitat)
• -6.47 (0.0015) odds of being in a higher category are 99% lower in south vs northern (for C habitat)

as_draws_df (posteriors)

hughes.brm3 |> as_draws_df()

# A draws_df: 800 iterations, 3 chains, and 835 variables
   b_Intercept[1] b_Intercept[2] b_Intercept[3] b_Intercept[4] b_HABITATF
1            -8.3           -5.2           -2.6          -0.35     -0.295
2            -8.7           -5.6           -2.9          -0.54     -0.339
3            -8.3           -5.2           -2.5          -0.24     -0.227
4            -8.4           -5.2           -2.5          -0.14      0.041
5            -8.5           -5.4           -2.8          -0.50     -0.440
6            -8.7           -5.6           -3.1          -0.63     -0.453
7            -8.7           -5.3           -2.7          -0.41     -0.440
8            -8.4           -5.3           -2.7          -0.22     -0.145
9            -8.3           -5.1           -2.7          -0.35     -0.339
10           -8.7           -5.5           -2.8          -0.29     -0.533
   b_HABITATL b_HABITATU b_SECTORCentral
1        -1.5       0.69           -0.93
2        -1.3       1.33           -1.76
3        -1.2       0.62           -1.44
4        -1.3       1.50           -2.30
5        -1.5       0.97           -2.01
6        -1.1       0.95           -2.34
7        -1.3       0.43           -1.59
8        -1.3       1.62           -1.49
9        -1.6       0.39           -1.77
10       -1.4       1.18           -2.00
# ... with 2390 more draws, and 827 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}

hughes.brm3 |>
  as_draws_df() |>
  exp() |>
  ## dplyr::select(matches("^b_Intercept.*|^b_HABITAT.*|^b_SECTOR.*")) |>
  dplyr::select(matches("^b_*")) |>
  summarise_draws(
    median,
    HDInterval::hdi,
    Pl = ~ mean(.x < 1),
    Pg = ~ mean(.x > 1),
    rhat,
    ess_bulk
  )

Warning: Dropping 'draws_df' class as required metadata was removed.

# A tibble: 15 × 8
   variable                   median   lower   upper    Pl     Pg  rhat ess_bulk
   <chr>                       <dbl>   <dbl>   <dbl> <dbl>  <dbl> <dbl>    <dbl>
 1 b_Intercept[1]           0.000208 1.22e-4 3.15e-4 1     0      1.000    2096.
 2 b_Intercept[2]           0.00477  3.12e-3 6.86e-3 1     0      1.00     2154.
 3 b_Intercept[3]           0.0647   4.51e-2 8.71e-2 1     0      1.00     2071.
 4 b_Intercept[4]           0.670    4.81e-1 8.95e-1 0.995 0.005  1.00     2108.
 5 b_HABITATF               0.713    4.86e-1 9.92e-1 0.970 0.0304 1.000    2254.
 6 b_HABITATL               0.261    1.81e-1 3.65e-1 1     0      1.00     2246.
 7 b_HABITATU               3.25     1.28e+0 5.96e+0 0     1      1.000    2383.
 8 b_SECTORCentral          0.179    8.33e-2 3.13e-1 1     0      1.00     1946.
 9 b_SECTORSouth            0.00158  8.75e-4 2.50e-3 1     0      1.00     1846.
10 b_HABITATF:SECTORCentral 0.524    2.50e-1 9.15e-1 0.98  0.02   1.00     2382.
11 b_HABITATL:SECTORCentral 1.23     6.13e-1 2.17e+0 0.244 0.756  1.000    2443.
12 b_HABITATU:SECTORCentral 0.558    1.22e-1 1.31e+0 0.879 0.121  1.00     2434.
13 b_HABITATF:SECTORSouth   0.668    3.88e-1 1.01e+0 0.958 0.0425 1.000    2172.
14 b_HABITATL:SECTORSouth   0.605    3.56e-1 9.16e-1 0.985 0.0154 1.00     2279.
15 b_HABITATU:SECTORSouth   0.700    2.62e-1 1.39e+0 0.826 0.174  1.000    2417.

\(R^2\)

hughes.brm3 |>
  bayes_R2(re.form = NA, summary = FALSE) |>
  median_hdci()

Warning: Predictions are treated as continuous variables in 'bayes_R2' which is
likely invalid for ordinal families.

          y      ymin      ymax .width .point .interval
1 0.6113859 0.5906381 0.6334466   0.95 median      hdci

hughes.brm3 |>
  bayes_R2(re.form = ~ (1 | REEF), summary = FALSE) |>
  median_hdci()

Warning: Predictions are treated as continuous variables in 'bayes_R2' which is
likely invalid for ordinal families.

          y      ymin      ymax .width .point .interval
1 0.8278381 0.8212367 0.8340643   0.95 median      hdci

12 Further investigations

brms

Cellmeans

Lets start by calculating the cell means

newdata <- with(hughes, expand.grid(
  HABITAT = levels(HABITAT),
  SECTOR = levels(SECTOR)
))
newdata

   HABITAT  SECTOR
1        C   North
2        F   North
3        L   North
4        U   North
5        C Central
6        F Central
7        L Central
8        U Central
9        C   South
10       F   South
11       L   South
12       U   South

predict(hughes.brm3, newdata = newdata, re_formula = NA)

          P(Y = 0)    P(Y = 1)   P(Y = 2)   P(Y = 3)    P(Y = 4)
 [1,] 0.0004166667 0.003750000 0.06291667 0.33000000 0.602916667
 [2,] 0.0004166667 0.006666667 0.07875000 0.41125000 0.502916667
 [3,] 0.0004166667 0.017083333 0.17583333 0.52166667 0.285000000
 [4,] 0.0000000000 0.001666667 0.02250000 0.15458333 0.821250000
 [5,] 0.0008333333 0.029166667 0.25416667 0.50083333 0.215000000
 [6,] 0.0016666667 0.083750000 0.41875000 0.39916667 0.096666667
 [7,] 0.0037500000 0.081666667 0.44041667 0.39208333 0.082083333
 [8,] 0.0004166667 0.019166667 0.15250000 0.48083333 0.347083333
 [9,] 0.1270833333 0.622500000 0.22791667 0.02000000 0.002500000
[10,] 0.2179166667 0.644583333 0.12208333 0.01375000 0.001666667
[11,] 0.4708333333 0.482916667 0.04250000 0.00250000 0.001250000
[12,] 0.0616666667 0.518333333 0.35666667 0.05791667 0.005416667

sum(0:4 * (c(0.0000000000, 0.0041666667, 0.04666667, 0.35000000, 0.5991666667)))

[1] 3.544167

add_epred_draws(hughes.brm3,
  newdata = newdata,
  re_formula = NA
) |>
  ## filter(.draw == 1)
  mutate(fit = as.numeric(as.character(.category)) * .epred) |>
  group_by(HABITAT, SECTOR, .draw) |>
  summarise(fit = sum(fit)) |>
  summarise_draws(
    median,
    HDInterval::hdi
  ) |>
  arrange(SECTOR, HABITAT)

`summarise()` has grouped output by 'HABITAT', 'SECTOR'. You can override using
the `.groups` argument.

# A tibble: 12 × 6
# Groups:   HABITAT, SECTOR [12]
   HABITAT SECTOR  variable median lower upper
   <fct>   <fct>   <chr>     <dbl> <dbl> <dbl>
 1 C       North   fit       3.53  3.45  3.63 
 2 F       North   fit       3.43  3.28  3.55 
 3 L       North   fit       3.06  2.91  3.23 
 4 U       North   fit       3.81  3.67  3.91 
 5 C       Central fit       2.92  2.69  3.13 
 6 F       Central fit       2.53  2.24  2.78 
 7 L       Central fit       2.47  2.18  2.74 
 8 U       Central fit       3.14  2.76  3.47 
 9 C       South   fit       1.16  1.05  1.27 
10 F       South   fit       0.933 0.804 1.05 
11 L       South   fit       0.599 0.482 0.724
12 U       South   fit       1.43  1.22  1.60 

hughes.cellmeans <-
  add_epred_draws(hughes.brm3, newdata = newdata, re_formula = NA) |>
  mutate(fit = as.numeric(as.character(.category)) * .epred) |>
  group_by(HABITAT, SECTOR, .draw) |>
  summarise(fit = sum(fit)) |>
  summarise_draws(
    median,
    HDInterval::hdi
  ) |>
  arrange(SECTOR, HABITAT)

`summarise()` has grouped output by 'HABITAT', 'SECTOR'. You can override using
the `.groups` argument.

fig1 <-
  hughes.cellmeans |>
  ggplot(aes(y = median, x = HABITAT)) +
  geom_hline(yintercept = 1, linetype = "dashed", size = 0.1) +
  geom_hline(yintercept = 2, linetype = "dashed", size = 0.1) +
  geom_hline(yintercept = 3, linetype = "dashed", size = 0.1) +
  geom_pointrange(aes(ymin = lower, ymax = upper)) +
  facet_grid(~SECTOR) +
  scale_y_continuous("Bleaching score", breaks = (0:4), labels = 0:4, limits = c(0, 4), expand = c(0, 0)) +
  theme_bw() +
  theme(panel.spacing.y = unit(10, "pt"))

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

## Alternative
## marginaleffects::avg_predictions(hughes.brm3, newdata = newdata, c("SECTOR", "HABITAT"), re_formula = NA) |>
##     marginaleffects::posterior_draws() |>
##   dplyr::select(SECTOR, HABITAT, group, estimate) |>
##     group_by(SECTOR, HABITAT, group) |>
##     summarise_draws(median)

## marginaleffects::avg_predictions(hughes.brm3, newdata = newdata, c("SECTOR", "HABITAT"), re_formula = NA) |>
##     marginaleffects::posterior_draws() |>
##     ggplot(aes(x = estimate, y = HABITAT, fill = group)) +
##   stat_halfeye() +
##   facet_grid(~SECTOR)
##     head()

Effects

newdata <- with(hughes, expand.grid(
  HABITAT = levels(HABITAT),
  SECTOR = levels(SECTOR)
))

hughes.eff <-
  add_epred_draws(hughes.brm3,
    newdata = newdata,
    re_formula = NA
  ) |>
  mutate(fit = as.numeric(as.character(.category)) * .epred) |>
  group_by(HABITAT, SECTOR, .draw) |>
  summarise(fit = log(sum(fit))) |>
  ## compare_levels(var = fit, by = HABITAT, comparison = "pairwise") |>
  tidybayes::compare_levels(
    var = fit, by = HABITAT,
    ## comparison = emmeans_comparison("tukey", reverse = TRUE)) |>
    comparison = emmeans_comparison("revpairwise")
    ## comparison = emmeans_comparison("pairwise")
    ## comparison = "pairwise"
  ) |>
  mutate(fit = exp(fit)) |>
  group_by(SECTOR, HABITAT) |>
  summarise_draws(
    median,
    HDInterval::hdi
  )

`summarise()` has grouped output by 'HABITAT', 'SECTOR'. You can override using
the `.groups` argument.

fig2 <-
  hughes.eff |>
  ggplot(aes(x = median, y = HABITAT)) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  geom_pointrange(aes(xmin = lower, xmax = upper)) +
  ## facet_grid(~SECTOR, scales = "free") +
  facet_grid(~SECTOR) +
  scale_x_continuous("Effect size (percentage change in bleaching category)",
    trans = scales::log2_trans(),
    breaks = seq(0.25, 2.5, by = 0.25),
    ## breaks = scales::trans_breaks("log10", function(x) exp(x)),
    ## breaks = scales::trans_breaks("log2", function(x) 2^x),
    ## labels = scales::trans_format("log2", function(x) x),
    ## labels = function(x) x
    labels = function(x) (x - 1) * 100
  ) +
  theme_bw()

fig1 + fig2

## Alternative with slab

hughes.eff <-
  add_epred_draws(hughes.brm3,
    newdata = newdata,
    re_formula = NA
  ) |>
  mutate(fit = as.numeric(as.character(.category)) * .epred) |>
  group_by(HABITAT, SECTOR, .draw) |>
  summarise(fit = log(sum(fit))) |>
  ## compare_levels(var = fit, by = HABITAT, comparison = "pairwise") |>
  tidybayes::compare_levels(
    var = fit, by = HABITAT,
    ## comparison = emmeans_comparison("tukey", reverse = TRUE)) |>
    comparison = emmeans_comparison("revpairwise")
    ## comparison = emmeans_comparison("pairwise")
    ## comparison = "pairwise"
  ) |>
  mutate(fit = exp(fit))

`summarise()` has grouped output by 'HABITAT', 'SECTOR'. You can override using
the `.groups` argument.

fig2 <-
  hughes.eff |>
  ggplot(aes(x = fit, y = HABITAT)) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  stat_halfeye(slab_fill = "orange", slab_alpha = 0.5, normalize = "panels") +
  facet_grid(~SECTOR, scales = "free") +
  ## facet_grid(~SECTOR) +
  scale_x_continuous("Effect size (percentage change in bleaching category)",
    trans = scales::log2_trans(),
    breaks = seq(0.25, 2.5, by = 0.25),
    ## breaks = scales::trans_breaks("log10", function(x) exp(x)),
    ## breaks = scales::trans_breaks("log2", function(x) 2^x),
    ## labels = scales::trans_format("log2", function(x) x),
    ## labels = function(x) x
    labels = function(x) (x - 1) * 100
  ) +
  theme_bw()

fig1 + fig2