Bayesian GLMM Part7

Author

Murray Logan

Published

09/09/2025

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(ggeffects)
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(tidyverse) # for data wrangling
library(brms)
library(tidybayes)
library(broom.mixed)
library(rstan)
library(cmdstanr)
library(patchwork)
library(DHARMa)
library(easystats)
library(modelsummary)
source("helperFunctions.R")

2 Scenario

In an honours thesis from (1992), Mullens was investigating the ways that cane toads ( Bufo marinus ) respond to conditions of hypoxia. Toads show two different kinds of breathing patterns, lung or buccal, requiring them to be treated separately in the experiment. Her aim was to expose toads to a range of O2 concentrations, and record their breathing patterns, including parameters such as the expired volume for individual breaths. It was desirable to have around 8 replicates to compare the responses of the two breathing types, and the complication is that animals are expensive, and different individuals are likely to have different O2 profiles (leading to possibly reduced power). There are two main design options for this experiment;

  • One animal per O2 treatment, 8 concentrations, 2 breathing types. With 8 replicates the experiment would require 128 animals, but that this could be analysed as a completely randomized design
  • One O2 profile per animal, so that each animal would be used 8 times and only 16 animals are required (8 lung and 8 buccal breathers)

Mullens decided to use the second option so as to reduce the number of animals required (on financial and ethical grounds). By selecting this option, she did not have a set of independent measurements for each oxygen concentration, by repeated measurements on each animal across the 8 oxygen concentrations.

Figure 1: Toad
Table 1: Format of mullens.csv data file
BREATH TOAD O2LEVEL FREQBUC SFREQBUC
lung a 0 10.6 3.256
lung a 5 18.8 4.336
lung a 10 17.4 4.171
lung a 15 16.6 4.074
... ... ... ... ...
Table 2: Description of the variables in the mullens data file
BREATH Categorical listing of the breathing type treatment (buccal = buccal breathing toads, lung = lung breathing toads). This is the between subjects (plots) effect and applies to the whole toads (since a single toad can only be one breathing type - either lung or buccal). Equivalent to Factor A (between plots effect) in a split-plot design
TOAD These are the subjects (equivalent to the plots in a split-plot design: Factor B). The letters in this variable represent the labels given to each individual toad.
O2LEVEL 0 through to 50 represent the the different oxygen concentrations (0% to 50%). The different oxygen concentrations are equivalent to the within plot effects in a split-plot (Factor C).
FREQBUC The frequency of buccal breathing - the response variable
SFREQBUC Square root transformed frequency of buccal breathing - the response variable

3 Read in the data

mullens <- read_csv("../data/mullens.csv", trim_ws = TRUE)
Rows: 168 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): BREATH, TOAD
dbl (3): O2LEVEL, FREQBUC, SFREQBUC

ℹ 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(mullens)
Rows: 168
Columns: 5
$ BREATH   <chr> "lung", "lung", "lung", "lung", "lung", "lung", "lung", "lung…
$ TOAD     <chr> "a", "a", "a", "a", "a", "a", "a", "a", "b", "b", "b", "b", "…
$ O2LEVEL  <dbl> 0, 5, 10, 15, 20, 30, 40, 50, 0, 5, 10, 15, 20, 30, 40, 50, 0…
$ FREQBUC  <dbl> 10.6, 18.8, 17.4, 16.6, 9.4, 11.4, 2.8, 4.4, 21.6, 17.4, 22.4…
$ SFREQBUC <dbl> 3.255764, 4.335897, 4.171331, 4.074310, 3.065942, 3.376389, 1…
## Explore the first 6 rows of the data
head(mullens)
# A tibble: 6 × 5
  BREATH TOAD  O2LEVEL FREQBUC SFREQBUC
  <chr>  <chr>   <dbl>   <dbl>    <dbl>
1 lung   a           0    10.6     3.26
2 lung   a           5    18.8     4.34
3 lung   a          10    17.4     4.17
4 lung   a          15    16.6     4.07
5 lung   a          20     9.4     3.07
6 lung   a          30    11.4     3.38
str(mullens)
spc_tbl_ [168 × 5] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ BREATH  : chr [1:168] "lung" "lung" "lung" "lung" ...
 $ TOAD    : chr [1:168] "a" "a" "a" "a" ...
 $ O2LEVEL : num [1:168] 0 5 10 15 20 30 40 50 0 5 ...
 $ FREQBUC : num [1:168] 10.6 18.8 17.4 16.6 9.4 11.4 2.8 4.4 21.6 17.4 ...
 $ SFREQBUC: num [1:168] 3.26 4.34 4.17 4.07 3.07 ...
 - attr(*, "spec")=
  .. cols(
  ..   BREATH = col_character(),
  ..   TOAD = col_character(),
  ..   O2LEVEL = col_double(),
  ..   FREQBUC = col_double(),
  ..   SFREQBUC = col_double()
  .. )
 - attr(*, "problems")=<externalptr> 
mullens |> datawizard::data_codebook()
mullens (168 rows and 5 variables, 5 shown)

ID | Name     | Type      | Missings |  Values |           N
---+----------+-----------+----------+---------+------------
1  | BREATH   | character | 0 (0.0%) |  buccal | 104 (61.9%)
   |          |           |          |    lung |  64 (38.1%)
---+----------+-----------+----------+---------+------------
2  | TOAD     | character | 0 (0.0%) |       a |   8 ( 4.8%)
   |          |           |          |       b |   8 ( 4.8%)
   |          |           |          |       c |   8 ( 4.8%)
   |          |           |          |       d |   8 ( 4.8%)
   |          |           |          |       e |   8 ( 4.8%)
   |          |           |          |       f |   8 ( 4.8%)
   |          |           |          |       g |   8 ( 4.8%)
   |          |           |          |       h |   8 ( 4.8%)
   |          |           |          |       i |   8 ( 4.8%)
   |          |           |          |       j |   8 ( 4.8%)
   |          |           |          |   (...) |            
---+----------+-----------+----------+---------+------------
3  | O2LEVEL  | numeric   | 0 (0.0%) | [0, 50] |         168
---+----------+-----------+----------+---------+------------
4  | FREQBUC  | numeric   | 0 (0.0%) | [0, 49] |         168
---+----------+-----------+----------+---------+------------
5  | SFREQBUC | numeric   | 0 (0.0%) |  [0, 7] |         168
------------------------------------------------------------
mullens |> modelsummary::datasummary_skim()
Unique Missing Pct. Mean SD Min Median Max Histogram
O2LEVEL 8 0 21.2 16.4 0.0 17.5 50.0
FREQBUC 90 0 13.7 10.6 0.0 10.1 49.0
SFREQBUC 90 0 3.4 1.5 0.0 3.2 7.0
N %
BREATH buccal 104 61.9
lung 64 38.1
TOAD a 8 4.8
b 8 4.8
c 8 4.8
d 8 4.8
e 8 4.8
f 8 4.8
g 8 4.8
h 8 4.8
i 8 4.8
j 8 4.8
k 8 4.8
l 8 4.8
m 8 4.8
n 8 4.8
o 8 4.8
p 8 4.8
q 8 4.8
r 8 4.8
s 8 4.8
t 8 4.8
u 8 4.8 
mullens |> modelsummary::datasummary_skim(by = c("BREATH", "O2LEVEL"))
BREATH O2LEVEL Unique Missing Pct. Mean SD Min Median Max Histogram
FREQBUC lung 0.0 6 0 4.3 5.2 0.0 2.7 13.8
lung 5.0 6 0 8.2 7.7 0.0 5.3 18.8
lung 10.0 8 0 13.2 7.7 0.0 16.7 22.4
lung 15.0 8 0 11.4 8.2 0.0 10.6 27.0
lung 20.0 8 0 11.9 8.3 4.4 9.6 31.0
lung 30.0 8 0 10.2 6.7 3.8 8.7 25.0
lung 40.0 7 0 11.1 15.5 2.8 6.7 49.0
lung 50.0 8 0 7.2 5.8 2.8 5.3 21.0
buccal 0.0 12 0 25.7 11.4 8.4 28.0 45.8
buccal 5.0 12 0 23.6 11.2 7.6 21.4 44.0
buccal 10.0 12 0 20.2 9.0 9.6 16.2 38.0
buccal 15.0 13 0 17.9 10.4 5.2 16.8 40.2
buccal 20.0 12 0 12.7 8.1 3.0 9.4 29.2
buccal 30.0 13 0 13.9 11.1 2.8 9.6 39.0
buccal 40.0 12 0 8.1 4.9 0.8 6.4 16.8
buccal 50.0 12 0 6.9 5.6 2.4 5.2 22.2
SFREQBUC lung 0.0 6 0 1.5 1.5 0.0 1.6 3.7
lung 5.0 6 0 2.4 1.7 0.0 2.3 4.3
lung 10.0 8 0 3.3 1.6 0.0 4.1 4.7
lung 15.0 8 0 3.1 1.5 0.0 3.2 5.2
lung 20.0 8 0 3.3 1.0 2.1 3.1 5.6
lung 30.0 8 0 3.1 1.0 1.9 2.9 5.0
lung 40.0 7 0 2.9 1.7 1.7 2.6 7.0
lung 50.0 8 0 2.6 0.9 1.7 2.3 4.6
buccal 0.0 12 0 4.9 1.2 2.9 5.3 6.8
buccal 5.0 12 0 4.7 1.2 2.8 4.6 6.6
buccal 10.0 12 0 4.4 1.0 3.1 4.0 6.2
buccal 15.0 13 0 4.1 1.2 2.3 4.1 6.3
buccal 20.0 12 0 3.4 1.1 1.7 3.1 5.4
buccal 30.0 13 0 3.5 1.4 1.7 3.1 6.2
buccal 40.0 12 0 2.7 0.9 0.9 2.5 4.1
buccal 50.0 12 0 2.5 0.9 1.5 2.3 4.7
N %
BREATH buccal 104 61.9
lung 64 38.1
TOAD a 8 4.8
b 8 4.8
c 8 4.8
d 8 4.8
e 8 4.8
f 8 4.8
g 8 4.8
h 8 4.8
i 8 4.8
j 8 4.8
k 8 4.8
l 8 4.8
m 8 4.8
n 8 4.8
o 8 4.8
p 8 4.8
q 8 4.8
r 8 4.8
s 8 4.8
t 8 4.8
u 8 4.8

4 Exploratory data analysis

Model formula: \[ \begin{align} y_i \sim{} \mathcal{Pois}(\lambda_i)\\ ln(\lambda_i) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \end{align} \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of copper, distance and their interaction on the number of number of worms. Area of the place segment was also incorporated as an offset. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual plates.

In this example, the individual TOADs are the random blocks. Each toad can only be of one BREATH type (they are either predominantly buccal breathers or predominantly lung breathers). Hence, BREATH is a between TOAD (block) effect. The frequency of buccal breathing of ach TOAD was measured under eight different oxygen levels and thus, these represent the within block effect, as will the interaction between breathing type and oxygen level.

The response in this case is a little tricky since it is a proportion, but without the full original measurements. For example a frequency of 50% would indicate that half of the breaths taken by the toad during the monitoring phase were buccal breaths. Ideally, it would be good to have the actual counts (both the number of buccal breaths and the total number of breaths). That way, we could model the data against a binomial distribution.

As it is, we only have the proportion (as a percentage). Although we could model this against a beta distribution, this could be complicated if there are proportions of either 0 or 1 (100).

To help guide us through this and the other typical model assumptions, lets start by graphing the frequency of buccal breathing against breathing type and oxygen level. However, before we do, we need to make sure that all categorical variables are declared as factors (including the random effect). If we intend to model against a beta distribution, we will need a version of the response that is represented on a scale between 0 and 1 (but not include 0 or 1).

mullens <- mullens |>
  mutate(
    BREATH = factor(BREATH),
    TOAD = factor(TOAD),
    pBUC = FREQBUC / 100,
    pzBUC = ifelse(pBUC == 0, 0.01, pBUC)
  )

So starting with the raw response.

ggplot(mullens, aes(y = FREQBUC, x = factor(O2LEVEL), color = BREATH)) +
  geom_boxplot()

Conclusions:

  • there is a very clear relationship between mean and variance (boxplots that are higher up the y-axis are taller).
  • it might be possible to address this via a logarithmic or root transformation, however this is the least favourable option (particularly root transformations due to issues with back-transformations).

If we intend to use a beta distribution, we could repeat the above with the scaled response. Since the rescaling maintains the same ranking and relative spacing, the plot will look the same as above, just with a different y-axis. However, our interpretation will change. Under a beta distribution (with logit link), we expect that the variance and mean will be related. We expect that the distributions will become more varied and more symmetrical as the expected mean shifts towards 0.5. Distributions approaching 0 and 1 will be more asymmetrical and smaller. This indeed does seem to be the case here.

Now lets explore the oxygen trends (separately for each breathing type).

ggplot(mullens, aes(y = pzBUC, x = O2LEVEL, color = BREATH)) +
  geom_smooth() +
  geom_point()
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Conclusions:

  • although linearity seems reasonable for the buccal breathers, there is definitely evidence of non linearity in the lung breathers.
  • it might be interesting to fit polynomial trends for oxygen concentration.
ggplot(mullens, aes(y = pzBUC, x = O2LEVEL, color = BREATH)) +
  geom_smooth() +
  geom_point() +
  facet_wrap(~ BREATH + TOAD, scales = "free")
`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

# facet_grid(TOAD~BREATH)

Conclusions:

  • it does appear that different toads (within the same breathing type) have very different levels of buccal breathing
  • the individual toads also have varying responses to changes in oxygen concentration
  • it might be useful to explore random intercept/slope models.

5 Fit the model

Explore the default priors

mullens.form <- bf(pBUC ~ BREATH * poly(O2LEVEL, 3) + (1 | TOAD),
  family = zero_one_inflated_beta()
)
get_prior(mullens.form, data = mullens)
                prior     class                     coef group resp dpar nlpar
               (flat)         b                                               
               (flat)         b               BREATHlung                      
               (flat)         b BREATHlung:polyO2LEVEL31                      
               (flat)         b BREATHlung:polyO2LEVEL32                      
               (flat)         b BREATHlung:polyO2LEVEL33                      
               (flat)         b            polyO2LEVEL31                      
               (flat)         b            polyO2LEVEL32                      
               (flat)         b            polyO2LEVEL33                      
           beta(1, 1)       coi                                               
 student_t(3, 0, 2.5) Intercept                                               
    gamma(0.01, 0.01)       phi                                               
 student_t(3, 0, 2.5)        sd                                               
 student_t(3, 0, 2.5)        sd                           TOAD                
 student_t(3, 0, 2.5)        sd                Intercept  TOAD                
           beta(1, 1)       zoi                                               
 lb ub       source
            default
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
  0  1      default
            default
  0         default
  0         default
  0    (vectorized)
  0    (vectorized)
  0  1      default

Start by fitting with priors only Remember that when fitting polynomials, they are essentially centered…

mullens |>
  group_by(BREATH) |>
  summarise(
    logit(median(pBUC)),
    logit(mad(pBUC))
  )
# A tibble: 2 × 3
  BREATH `logit(median(pBUC))` `logit(mad(pBUC))`
  <fct>                  <dbl>              <dbl>
1 buccal                 -1.87              -2.05
2 lung                   -2.48              -2.79
median(logit(mullens$pBUC))
Warning in logit(mullens$pBUC): proportions remapped to (0.025, 0.975)
[1] -1.983737
mad(logit(mullens$pBUC))
Warning in logit(mullens$pBUC): proportions remapped to (0.025, 0.975)
[1] 0.8044649
sd(logit(mullens$pBUC)) / apply(model.matrix(~ BREATH * poly(O2LEVEL, 3), data = mullens), 2, sd)
Warning in logit(mullens$pBUC): proportions remapped to (0.025, 0.975)
                 (Intercept)                   BREATHlung 
                         Inf                     1.641316 
           poly(O2LEVEL, 3)1            poly(O2LEVEL, 3)2 
                   10.331042                    10.331042 
           poly(O2LEVEL, 3)3 BREATHlung:poly(O2LEVEL, 3)1 
                   10.331042                    16.738202 
BREATHlung:poly(O2LEVEL, 3)2 BREATHlung:poly(O2LEVEL, 3)3 
                   16.738202                    16.738202 
standist::visualize("gamma(0.01, 0.01)", "gamma(2,1)", "gamma(1,1)", xlim = c(0, 10))

standist::visualize("beta(1,1)", xlim = c(0, 1))

priors <- prior(normal(-2, 2), class = "Intercept") +
  prior(normal(0, 1), class = "b") +
  ## prior(normal(0,1), class='b', coef = "BREATHlung") +
  # prior(gamma(2,1), class='sd') +
  prior(student_t(3, 0, 2), class = "sd") +
  prior(gamma(0.01, 0.01), class = "phi") +
  ## prior(gamma(2, 1), class='phi') +
  prior(beta(1, 1), class = "zoi") +
  prior(beta(1, 1), class = "coi")

mullens.form <- bf(pBUC ~ BREATH * poly(O2LEVEL, 3) + (1 | TOAD),
  family = zero_one_inflated_beta()
)

mullens.brm <- brm(mullens.form,
  data = mullens,
  prior = priors,
  sample_prior = "only",
  iter = 5000,
  warmup = 2500,
  thin = 5,
  chains = 3, cores = 3,
  refresh = 0,
  seed = 123,
  control = list(adapt_delta = 0.99, max_treedepth = 20),
  backend = "rstan"
)
Compiling Stan program...
Start sampling
Warning: There were 31 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.
Warning: Examine the pairs() plot to diagnose sampling problems
mullens.brm |>
  conditional_effects(effects = "O2LEVEL:BREATH") |>
  plot(points = TRUE)

mullens.brm2 <- update(mullens.brm,
  sample_prior = "yes", refresh = 0,
  cores = 3, seed = 123
)
The desired updates require recompiling the model
Compiling Stan program...
Start sampling
priors <- prior(normal(-2, 2), class = "Intercept") +
  prior(normal(0, 2), class = "b") +
  prior(student_t(3, 0, 2), class = "sd") +
  prior(gamma(0.01, 0.01), class = "phi") +
  prior(beta(1, 1), class = "zoi") +
  prior(beta(1, 1), class = "coi")
mullens.brm3 <- update(mullens.brm2,
  prior = priors, refresh = 0,
  cores = 3, seed = 123
)
The desired updates require recompiling the model
Compiling Stan program...
Start sampling
mullens.brm2 |>
  conditional_effects(effects = "O2LEVEL:BREATH") |>
  plot(points = TRUE)

mullens.brm2 |>
  ggpredict(~ O2LEVEL | BREATH) |>
  plot(show_data = TRUE, jitter = FALSE)
You are calculating adjusted predictions on the population-level (i.e.
  `type = "fixed"`) for a *generalized* linear mixed model.
  This may produce biased estimates due to Jensen's inequality. Consider
  setting `bias_correction = TRUE` to correct for this bias.
  See also the documentation of the `bias_correction` argument.
Note: uncertainty of error terms are not taken into account. Consider
  setting `interval` to "prediction". This will call `posterior_predict()`
  instead of `posterior_epred()`.
Warning in glmmTMB::glmmTMB(formula = cond_formula, family = family, data =
data): some components missing from 'family': downstream methods may fail
Warning in glmmTMB::glmmTMB(formula = cond_formula, family = family, data =
data): some components missing from 'family': downstream methods may fail
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

mullens.brm3 |>
  conditional_effects(effects = "O2LEVEL:BREATH") |>
  plot(points = TRUE)

mullens.brm2 |> SUYR_prior_and_posterior()

mullens.brm3 |> SUYR_prior_and_posterior()

MCMC sampling diagnostics ::: {.panel-tabset} ## brms :::: {.panel-tabset} ### stan plots

mullens.brm3$fit |> stan_trace()
mullens.brm3$fit |> stan_ac()
mullens.brm3$fit |> stan_rhat()
mullens.brm3$fit |> stan_ess()

:::: :::

6 Model validation

mullens.brm3 |> pp_check(type = "dens_overlay", ndraws = 200)

mullens.brm3 |> pp_check(group = "BREATH", type = "intervals_grouped")
mullens.brm3 |> pp_check(group = "O2LEVEL", type = "intervals_grouped")
mullens.resids <- make_brms_dharma_res(mullens.brm3, integerResponse = FALSE)
wrap_elements(~ testUniformity(mullens.resids)) +
  wrap_elements(~ plotResiduals(mullens.resids, form = factor(rep(1, nrow(mullens))))) +
  wrap_elements(~ plotResiduals(mullens.resids, quantreg = TRUE)) +
  wrap_elements(~ testDispersion(mullens.resids))

7 Partial effects plot

mullens.brm3 |>
  conditional_effects(effects = "O2LEVEL:BREATH") |>
  plot(points = TRUE)

7.1 ggpredict

mullens.brm3 |>
  ggpredict(~ O2LEVEL | BREATH) |>
  plot(show_data = TRUE, jitter = FALSE)
Note: uncertainty of error terms are not taken into account. Consider
  setting `interval` to "prediction". This will call `posterior_predict()`
  instead of `posterior_epred()`.
Warning in glmmTMB::glmmTMB(formula = cond_formula, family = family, data =
data): some components missing from 'family': downstream methods may fail
Warning in glmmTMB::glmmTMB(formula = cond_formula, family = family, data =
data): some components missing from 'family': downstream methods may fail
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

7.2 ggemmeans

mullens.brm3 |>
  ggemmeans(~ O2LEVEL | BREATH) |>
  plot(show_data = TRUE, jitter = FALSE)
Warning in glmmTMB::glmmTMB(formula = cond_formula, family = family, data =
data): some components missing from 'family': downstream methods may fail
Warning in glmmTMB::glmmTMB(formula = cond_formula, family = family, data =
data): some components missing from 'family': downstream methods may fail
Warning in glmmTMB::glmmTMB(formula = cond_formula, family = family, data =
data): some components missing from 'family': downstream methods may fail
Warning in glmmTMB::glmmTMB(formula = cond_formula, family = family, data =
data): some components missing from 'family': downstream methods may fail
Data points may overlap. Use the `jitter` argument to add some amount of
  random variation to the location of data points and avoid overplotting.

newdata <- with(mullens, list(
  O2LEVEL = seq(min(O2LEVEL), max(O2LEVEL), len = 100),
  BREATH = levels(BREATH)
))

mullens.em <-
  mullens.brm3 |>
  emmeans(~ O2LEVEL | BREATH, at = newdata, type = "response") |>
  as.data.frame()
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
mullens.em |>
  ggplot(aes(y = response, x = O2LEVEL, colour = BREATH, fill = BREATH)) +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD),
    alpha = 0.2, colour = NA
  ) +
  geom_line()

8 Model investigation

mullens.brm2 |> summary()
 Family: zero_one_inflated_beta 
  Links: mu = logit; phi = identity; zoi = identity; coi = identity 
Formula: pBUC ~ BREATH * poly(O2LEVEL, 3) + (1 | TOAD) 
   Data: mullens (Number of observations: 168) 
  Draws: 3 chains, each with iter = 5000; warmup = 2500; thin = 5;
         total post-warmup draws = 1500

Multilevel Hyperparameters:
~TOAD (Number of levels: 21) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.54      0.11     0.38     0.78 1.00     1226     1349

Regression Coefficients:
                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                   -1.77      0.16    -2.09    -1.46 1.00     1121
BREATHlung                  -0.38      0.25    -0.86     0.12 1.00     1303
polyO2LEVEL31               -4.34      0.50    -5.29    -3.33 1.00     1549
polyO2LEVEL32                0.25      0.48    -0.68     1.16 1.00     1438
polyO2LEVEL33                0.10      0.47    -0.80     1.01 1.00     1543
BREATHlung:polyO2LEVEL31     1.62      0.77     0.13     3.04 1.00     1401
BREATHlung:polyO2LEVEL32    -1.28      0.74    -2.72     0.15 1.00     1511
BREATHlung:polyO2LEVEL33     0.42      0.72    -0.95     1.84 1.00     1435
                         Tail_ESS
Intercept                    1141
BREATHlung                   1457
polyO2LEVEL31                1496
polyO2LEVEL32                1460
polyO2LEVEL33                1413
BREATHlung:polyO2LEVEL31     1486
BREATHlung:polyO2LEVEL32     1383
BREATHlung:polyO2LEVEL33     1497

Further Distributional Parameters:
    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi    34.18      4.44    25.97    43.57 1.00     1466     1333
zoi     0.05      0.02     0.02     0.08 1.00     1530     1313
coi     0.11      0.10     0.00     0.35 1.00     1407     1314

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).
mullens.brm3 |>
  as_draws_df() |>
  mutate(across(everything(), exp)) |>
  dplyr::select(matches("^b_.*|^sd_.*")) |>
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat,
    length,
    ess_bulk, ess_tail,
    Pl = ~ mean(.x < 1),
    Pg = ~ mean(.x > 1)
  ) |>
  knitr::kable()
Warning: Dropping 'draws_df' class as required metadata was removed.
variable median lower upper rhat length ess_bulk ess_tail Pl Pg
b_Intercept 0.1670042 0.1189784 0.2250192 1.0018252 1500 962.8235 1102.214 1.0000000 0.0000000
b_BREATHlung 0.6836448 0.3720328 1.0847048 1.0003004 1500 1051.6154 1110.226 0.9366667 0.0633333
b_polyO2LEVEL31 0.0031985 0.0007305 0.0084637 0.9997699 1500 1290.7153 1530.815 1.0000000 0.0000000
b_polyO2LEVEL32 1.5942031 0.3869634 3.8228459 0.9999437 1500 1503.3559 1280.052 0.1880000 0.8120000
b_polyO2LEVEL33 1.0841982 0.2266756 2.5565466 1.0004433 1500 1462.1332 1307.588 0.4380000 0.5620000
b_BREATHlung:polyO2LEVEL31 46.1674803 0.9903850 226.9142129 0.9996964 1500 1528.7125 1340.766 0.0006667 0.9993333
b_BREATHlung:polyO2LEVEL32 0.1092211 0.0071674 0.5038358 1.0070014 1500 1677.8320 1403.666 0.9953333 0.0046667
b_BREATHlung:polyO2LEVEL33 2.1747581 0.1256471 11.6626643 1.0002384 1500 1515.7266 1524.102 0.2246667 0.7753333
sd_TOAD__Intercept 1.7069957 1.4309808 2.1585563 1.0001341 1500 948.3748 1270.547 0.0000000 1.0000000 
mean(mullens$O2LEVEL)
[1] 21.25
0.163 / (1 + 0.163)
[1] 0.1401548
mullens.brm3 |>
  as_draws_df() |>
  dplyr::select(matches("^b_.*")) |>
  exp() |>
  summarise_draws(
    median,
    HDInterval::hdi,
    rhat,
    length,
    ess_bulk, ess_tail,
    Pl = ~ mean(.x < 1),
    Pg = ~ mean(.x > 1)
  )
Warning: Dropping 'draws_df' class as required metadata was removed.
# A tibble: 8 × 10
  variable         median   lower   upper  rhat length ess_bulk ess_tail      Pl
  <chr>             <dbl>   <dbl>   <dbl> <dbl>  <dbl>    <dbl>    <dbl>   <dbl>
1 b_Intercept     1.67e-1 1.19e-1 2.25e-1 1.00    1500     963.    1102. 1   e+0
2 b_BREATHlung    6.84e-1 3.72e-1 1.08e+0 1.00    1500    1052.    1110. 9.37e-1
3 b_polyO2LEVEL31 3.20e-3 7.30e-4 8.46e-3 1.000   1500    1291.    1531. 1   e+0
4 b_polyO2LEVEL32 1.59e+0 3.87e-1 3.82e+0 1.000   1500    1503.    1280. 1.88e-1
5 b_polyO2LEVEL33 1.08e+0 2.27e-1 2.56e+0 1.00    1500    1462.    1308. 4.38e-1
6 b_BREATHlung:p… 4.62e+1 9.90e-1 2.27e+2 1.000   1500    1529.    1341. 6.67e-4
7 b_BREATHlung:p… 1.09e-1 7.17e-3 5.04e-1 1.01    1500    1678.    1404. 9.95e-1
8 b_BREATHlung:p… 2.17e+0 1.26e-1 1.17e+1 1.00    1500    1516.    1524. 2.25e-1
# ℹ 1 more variable: Pg <dbl>
mullens.brm2 |> performance::r2()
# Bayesian R2 with Compatibility Interval

  Conditional R2: 0.553 (95% CI [0.452, 0.632])
     Marginal R2: 0.154 (95% CI [0.062, 0.287])
mullens.brm2 |>
  performance::r2_posterior() |>
  as.data.frame() |>
  median_hdci()
   R2_Bayes R2_Bayes.lower R2_Bayes.upper R2_Bayes_marginal
1 0.5530387      0.4617761      0.6421551         0.1535782
  R2_Bayes_marginal.lower R2_Bayes_marginal.upper .width .point .interval
1              0.06248566               0.2872945   0.95 median      hdci

9 Further investigations

mullens.brm3 |> emtrends(specs = "BREATH", var = "O2LEVEL", max.degree = 3)
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
degree = linear:
 BREATH O2LEVEL.trend lower.HPD upper.HPD
 buccal     -2.91e-02 -4.13e-02 -1.61e-02
 lung       -1.41e-02 -3.50e-02  4.50e-03

degree = quadratic:
 BREATH O2LEVEL.trend lower.HPD upper.HPD
 buccal      1.44e-04 -2.94e-04  5.07e-04
 lung       -7.83e-04 -1.48e-03 -1.06e-04

degree = cubic:
 BREATH O2LEVEL.trend lower.HPD upper.HPD
 buccal      1.99e-06 -2.58e-05  2.51e-05
 lung        2.08e-05 -2.26e-05  6.09e-05

Point estimate displayed: median 
HPD interval probability: 0.95 
mullens.brm3 |>
  emtrends(specs = "BREATH", var = "O2LEVEL", max.degree = 3) |>
  gather_emmeans_draws() |>
  summarise(median_hdci(.value),
    Pl = mean(.value < 0),
    Pg = mean(.value > 0)
  )
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
`summarise()` has grouped output by 'BREATH'. You can override using the
`.groups` argument.
# A tibble: 6 × 10
# Groups:   BREATH [2]
  BREATH degree         y     ymin     ymax .width .point .interval    Pl     Pg
  <fct>  <fct>      <dbl>    <dbl>    <dbl>  <dbl> <chr>  <chr>     <dbl>  <dbl>
1 buccal linear  -2.91e-2 -4.13e-2 -1.61e-2   0.95 median hdci      1     0     
2 buccal quadra…  1.44e-4 -2.87e-4  5.17e-4   0.95 median hdci      0.254 0.746 
3 buccal cubic    1.99e-6 -2.20e-5  2.89e-5   0.95 median hdci      0.438 0.562 
4 lung   linear  -1.41e-2 -3.54e-2  4.15e-3   0.95 median hdci      0.913 0.0867
5 lung   quadra… -7.83e-4 -1.50e-3 -1.18e-4   0.95 median hdci      0.985 0.0147
6 lung   cubic    2.08e-5 -2.26e-5  6.09e-5   0.95 median hdci      0.167 0.833 
mullens.brm3 |>
  emmeans(~ BREATH | O2LEVEL, at = list(O2LEVEL = c(0, 21, 50))) |>
  pairs() |>
  gather_emmeans_draws() |>
  mutate(.value = exp(.value)) |>
  summarise(median_hdci(.value),
    Pl = mean(.value < 1),
    Pg = mean(.value > 1)
  )
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
`summarise()` has grouped output by 'contrast'. You can override using the
`.groups` argument.
# A tibble: 3 × 10
# Groups:   contrast [1]
  contrast      O2LEVEL     y  ymin  ymax .width .point .interval       Pl    Pg
  <fct>           <dbl> <dbl> <dbl> <dbl>  <dbl> <chr>  <chr>        <dbl> <dbl>
1 buccal - lung       0  2.97 1.36   5.19   0.95 median hdci      0.000667 0.999
2 buccal - lung      21  1.19 0.604  1.92   0.95 median hdci      0.257    0.743
3 buccal - lung      50  1.06 0.522  1.91   0.95 median hdci      0.424    0.576
mullens.grid <- with(
  mullens,
  list(
    O2LEVEL = seq(min(O2LEVEL), max(O2LEVEL), len = 1000),
    BREATH = levels(BREATH)
  )
)
## mullens.grid <- with(mullens,
##   list(O2LEVEL = modelr::seq_range(O2LEVEL, n = 1000),
##     BREATH = levels(BREATH)))

## At what point does the evidence for the curves being different stop
mullens.brm3 |>
  emmeans(~ BREATH | O2LEVEL, at = mullens.grid) |>
  pairs() |>
  gather_emmeans_draws() |>
  mutate(.value = exp(.value)) |>
  summarise(median_hdci(.value),
    Pl = mean(.value < 1),
    Pg = mean(.value > 1)
  ) |>
  filter(Pg < 0.90) |>
  slice(1:3)
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
`summarise()` has grouped output by 'contrast'. You can override using the
`.groups` argument.
# A tibble: 3 × 10
# Groups:   contrast [1]
  contrast      O2LEVEL     y  ymin  ymax .width .point .interval    Pl    Pg
  <fct>           <dbl> <dbl> <dbl> <dbl>  <dbl> <chr>  <chr>     <dbl> <dbl>
1 buccal - lung    13.6  1.43 0.757  2.31   0.95 median hdci      0.101 0.899
2 buccal - lung    13.6  1.42 0.756  2.30   0.95 median hdci      0.101 0.899
3 buccal - lung    13.7  1.42 0.764  2.31   0.95 median hdci      0.103 0.897
## Find the peak
newdata <- mullens.brm3 |>
  emmeans(~ O2LEVEL | BREATH, at = mullens.grid) |>
  as.data.frame()
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
newdata |>
  group_by(BREATH) |>
  summarise(value = O2LEVEL[which.max(emmean)])
# A tibble: 2 × 2
  BREATH value
  <fct>  <dbl>
1 buccal   0  
2 lung    14.1
## With uncertainty
mullens.brm3 |>
  emmeans(~ O2LEVEL | BREATH, at = mullens.grid) |>
  gather_emmeans_draws() |>
  group_by(.draw, BREATH) |>
  summarise(value = O2LEVEL[which.max(.value)]) |>
  ungroup() |>
  group_by(BREATH) |>
  summarise(median_hdci(value))
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
`summarise()` has grouped output by '.draw'. You can override using the
`.groups` argument.
# A tibble: 2 × 7
  BREATH     y  ymin  ymax .width .point .interval
  <fct>  <dbl> <dbl> <dbl>  <dbl> <chr>  <chr>    
1 buccal   0       0   0     0.95 median hdci     
2 lung    14.1     0  22.8   0.95 median hdci     
# Or via Derivatives
## mullens.brm3 |> estimate_slopes(trend =  "O2LEVEL",
##   at = c("BREATH='lung'","O2LEVEL"))

## mullens.brm2 |>
##   emtrends(~O2LEVEL|BREATH, var = "O2LEVEL",
##     cov.red = \(x) seq(min(x), max(x), length = 10)
##   ) |>
##   summary() |>
##   ggplot(aes(y = O2LEVEL.trend, x = O2LEVEL, colour =  BREATH)) +
##   geom_hline(yintercept = 0, linetype = "dashed") +
##   geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = BREATH), alpha = 0.2) +
##   geom_line()

10 Summary figure

mullens.grid <- with(
  mullens,
  list(
    BREATH = levels(BREATH),
    O2LEVEL = modelr::seq_range(O2LEVEL, n = 100)
  )
)
newdata <- mullens.brm2 |>
  emmeans(~ O2LEVEL | BREATH, at = mullens.grid, type = "response") |>
  as.data.frame()
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
head(newdata)
   O2LEVEL BREATH  response lower.HPD upper.HPD
 0.0000000 buccal 0.2107378 0.1527691 0.2672369
 0.5050505 buccal 0.2091359 0.1548155 0.2668167
 1.0101010 buccal 0.2074718 0.1509585 0.2608619
 1.5151515 buccal 0.2055807 0.1562472 0.2643178
 2.0202020 buccal 0.2037584 0.1542856 0.2605384
 2.5252525 buccal 0.2019476 0.1535954 0.2582075

Point estimate displayed: median 
Results are back-transformed from the logit scale 
HPD interval probability: 0.95 
ggplot() +
  geom_ribbon(
    data = newdata,
    aes(
      ymin = lower.HPD, ymax = upper.HPD,
      x = O2LEVEL, fill = BREATH
    ), alpha = 0.3
  ) +
  geom_line(
    data = newdata,
    aes(y = response, x = O2LEVEL, color = BREATH)
  ) +
  scale_y_continuous("Buccal breathing rate", labels = function(x) 100 * x) +
  theme_classic()

## prior_summary(mullens.brm)


## pars <- mullens.brm |> get_variables()
## wch <- grepl('^b_.*|^sd_.*|phi', pars, perl=TRUE)

## g <- vector('list', length=sum(wch)-1)
## names(g) <- pars[wch][-1]
## for (i in pars[wch]) {
##     print(i)
##     if (i == 'b_Intercept') next
##     p <- mullens.brm |> hypothesis(paste0(i,'=0'), class='') |> plot()
##     g[[i]] <- p[[1]]
## }
## patchwork::wrap_plots(g)

stan_trace(mullens.brm$fit, pars = pars[wch])
stan_ac(mullens.brm$fit, pars = pars[wch])
stan_rhat(mullens.brm$fit, pars = pars[wch])
stan_rhat(mullens.brm$fit)
stan_ess(mullens.brm$fit)


preds <- posterior_predict(mullens.brm, nsamples = 250, summary = FALSE)
mullens.resids <- createDHARMa(
  simulatedResponse = t(preds),
  observedResponse = mullens$pBUC,
  fittedPredictedResponse = apply(preds, 2, median),
  integerResponse = FALSE
)
plot(mullens.resids)
testDispersion(mullens.resids)

11 References