library(tidyverse) # for data wrangling
library(car) # for regression diagnostics
library(broom) # for tidy output
library(ggfortify) # for model diagnostics
library(knitr) # for kable
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(brms)
library(broom.mixed)
library(tidybayes)
library(bayesplot)
library(standist) # for visualizing distributions
library(rstanarm)
library(cmdstanr) # for cmdstan
library(ggeffects)
library(rstan)
library(DHARMa)
library(ggridges)
library(easystats) # framework for stats, modelling and visualisation
library(patchwork)
library(modelsummary)
source("helperFunctions.R")Bayesian GLMM Part4
1 Preparations
Load the necessary libraries
2 Scenario
To investigate synergistic coral defence by mutualist crustaceans, Mckeon et al. (2012) conducted an aquaria experiment in which colonies of a coral species were placed in a tank along with a predatory sea star and one of four symbiont combinations:
- no symbiont,
- a crab symbiont
- a shrimp symbiont
- both a crab and shrimp symbiont.
The experiments were conducted in a large octagonal flow-through seawater tank that was partitioned into eight sections, which thereby permitted two of each of the four symbiont combinations to be observed concurrently. The tank was left overnight and in the morning, the presence of feeding scars on each coral colony was scored as evidence of predation. The experiments were repeated ten times, each time with fresh coral colonies, sea stars and symbiont.
The ten experimental times represent blocks (random effects) within which the symbiont type (fixed effect) are nested.
3 Read in the data
Rows: 80 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): SYMBIONT
dbl (2): BLOCK, PREDATION
ℹ 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.spc_tbl_ [80 × 3] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
$ BLOCK : num [1:80] 1 1 2 2 3 3 4 4 5 5 ...
$ PREDATION: num [1:80] 0 1 1 1 1 1 1 1 1 1 ...
$ SYMBIONT : chr [1:80] "none" "none" "none" "none" ...
- attr(*, "spec")=
.. cols(
.. BLOCK = col_double(),
.. PREDATION = col_double(),
.. SYMBIONT = col_character()
.. )
- attr(*, "problems")=<externalptr> mckeon (80 rows and 3 variables, 3 shown)
ID | Name | Type | Missings | Values | N
---+-----------+-----------+----------+---------+-----------
1 | BLOCK | numeric | 0 (0.0%) | [1, 10] | 80
---+-----------+-----------+----------+---------+-----------
2 | PREDATION | numeric | 0 (0.0%) | 0 | 30 (37.5%)
| | | | 1 | 50 (62.5%)
---+-----------+-----------+----------+---------+-----------
3 | SYMBIONT | character | 0 (0.0%) | both | 20 (25.0%)
| | | | crabs | 20 (25.0%)
| | | | none | 20 (25.0%)
| | | | shrimp | 20 (25.0%)
------------------------------------------------------------Since the response here is the presence or absence of predation (feeding scars), a binomial distribution is appropriate.
We need to make sure that the categorical variable and the random effect are defined as factors. When doing so, it might be valuable to rearrange the order of the fixed effect (SYMBIONT) such that the none group is considered the first group. This way, the other levels will all naturally be compared to this level (hence it will be treated as a reference of control group).
4 Exploratory data analysis
Random intercepts model formula: \[ \begin{align} y_i &\sim{} \mathcal{N}(n, p_i)\\ ln\left(\frac{p_i}{1-p_1}\right) &=\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i}\\ \beta_0 &\sim{} \mathcal{N}(0, 2.5)\\ \beta_1 &\sim{} \mathcal{N}(0, 3)\\ \boldsymbol{\gamma_j} &\sim{} \mathcal{N}(0, \sigma^2_b)\\ \sigma^2_b &\sim{} t(3, 0, 1.5)\\ \end{align} \]
Random intercept/slope model formula: \[ \begin{align} y_i &\sim{} \mathcal{N}(n, p_i)\\ ln\left(\frac{p_i}{1-p_1}\right) &=\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i}\\ \beta_0 &\sim{} \mathcal{N}(0, 2)\\ \beta_1 &\sim{} \mathcal{N}(0, 10)\\ \boldsymbol{\gamma_j} &\sim{} \mathcal{N}(0, \boldsymbol{\Sigma})\\ \boldsymbol{\Sigma} &= \boldsymbol{D}({\sigma_b^2})\boldsymbol{\Omega}\boldsymbol{D}({\sigma_b^2})\\ \boldsymbol{\Omega} &\sim{} LKJ(\zeta)\\ \sigma_b^2 &\sim{} \mathcal{t}(3,0,1.5)\\ \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 symbionts on the probability of the colony experiencing predation. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual coral colonies.
5 Fit the model
In rstanarm, 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.
Priors for model '.'
------
Intercept (after predictors centered)
~ normal(location = 0, scale = 2.5)
Coefficients
Specified prior:
~ normal(location = [0,0,0], scale = [2.5,2.5,2.5])
Adjusted prior:
~ normal(location = [0,0,0], scale = [5.74,5.74,5.74])
Covariance
~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
------
See help('prior_summary.stanreg') for more detailsThis tells us:
for the intercept (when the family is binomial), a normal prior with a mean of 0 and a standard deviation of 2.5 is used. These are the default priors for bernoulli models. A mean of 0 on a logit scale corresponds to a probability of 0.5. Hence the prior expected value is 0.5.
for the coefficients (in this case, just the slope), the default prior is a normal prior centred around 0 with a standard deviation of 2.5. For binomial models, this is then adjusted for the scale of the data by dividing the 2.5 by the standard deviation of the predictor (then rounded).
Lets now run with priors only so that we can explore the range of values they allow in the posteriors.
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()`.
Conclusions:
- it is very difficult to assess the priors from the predicted posteriors when the posteriors are on the response scale as they will always approach values of 0 and 1 (no matter how wide they are on the link scale).
Interestingly, if we instead use ggemmeans, it will (perhaps erroneously) generate the partial effects on the link scale (yet with an incorrectly labelled y-axis). Probability scale values of 0.99 and 0.01 (very large and small respectively) correspond to value of approximately 4.6 and -4.6 respectively on the logit scale.
In the following partial plot, the routine attempts to convert the y-axes tick marks into percentages. Hence a value of 0.1 is converted to 10%. Consequently, values of -5 and 5 are labelled as -500% and 500% respectively. We know that on the probability scale, values must asymptote towards 0 and 1. Posterior intervals beyond -5 and 5 might be considered wide.
Alternatively, we could create the plot ourselves…
mckeon_rstanarm1 %>%
emmeans(~SYMBIONT, type = "link") %>%
as.data.frame() %>%
ggplot(aes(y = emmean, x = SYMBIONT)) +
geom_hline(yintercept = c(5, -5), linetype = "dashed") +
geom_pointrange(aes(ymin = lower.HPD, ymax = upper.HPD)) +
geom_point(
data = mckeon, aes(y = PREDATION),
position = position_jitter(width = 0.2, height = 0),
alpha = 0.4, color = "red"
)
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 0 (since this is 0.5 on the logit scale) with a standard deviation of 2
- mean of 0: since
logit(0.5) - alternatively, an argument could be made for a mean of 0.51: since
logit(mean(mckeon$PREDATION)) - sd of 2: since this not too large on logit scale
- mean of 0: since
- \(\beta_1\): normal centred at 0 with a standard deviation of 0.5
- sd of 0.5: since
model.matrix(~SYMBIONT, data=mckeon)[,-1] %>% apply(2,sd)
- sd of 0.5: since
- \(\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.
I will also overlay the raw data for comparison.
mckeon_rstanarm2 <- stan_glmer(PREDATION ~ SYMBIONT + (1 | BLOCK),
data = mckeon,
family = binomial(link = "logit"),
prior_intercept = normal(0, 2.5, autoscale = FALSE),
prior = normal(0, 6, autoscale = FALSE),
prior_covariance = decov(1, 1, 1, 1),
prior_PD = TRUE,
iter = 5000,
warmup = 1000,
chains = 3,
thin = 5,
refresh = 0
)Note: uncertainty of error terms are not taken into account. Consider
setting `interval` to "prediction". This will call `posterior_predict()`
instead of `posterior_epred()`.
Now lets refit, conditioning on the data.
posterior_vs_prior(mckeon_rstanarm3,
color_by = "vs", group_by = TRUE,
facet_args = list(scales = "free_y")
)
Drawing from prior...
Conclusions:
- in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
Note: uncertainty of error terms are not taken into account. Consider
setting `interval` to "prediction". This will call `posterior_predict()`
instead of `posterior_epred()`.
Data points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
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.
mckeon_form <- bf(PREDATION | trials(1) ~ SYMBIONT + (1|BLOCK),
family=binomial(link='logit'))
#OR
mckeon_form <- bf(PREDATION ~ SYMBIONT + (1|BLOCK),
family=bernoulli(link='logit'))
options(width=150)
mckeon_form |> get_prior(data = mckeon) prior class coef group resp dpar nlpar lb ub source
(flat) b default
(flat) b SYMBIONTboth (vectorized)
(flat) b SYMBIONTcrabs (vectorized)
(flat) b SYMBIONTshrimp (vectorized)
student_t(3, 0, 2.5) Intercept default
student_t(3, 0, 2.5) sd 0 default
student_t(3, 0, 2.5) sd BLOCK 0 (vectorized)
student_t(3, 0, 2.5) sd Intercept BLOCK 0 (vectorized)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 0 (since this is 0.5 on the logit scale) with a standard deviation of 2.5
- mean of 0: since
logit(0.5) - sd of 2: since this not too large on logit scale
- alternatively, we could potentially use
logistic(0,1)
- mean of 0: since
- \(\beta_1\): normal centred at 0 with a standard deviation of 6
- sd of 6: since
2.5/model.matrix(~SYMBIONT, data=mckeon) %>% apply(2,sd)
- sd of 6: since
- \(\sigma_j\): half-cauchy with parameters 0 and 2.
- \(\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.
mckeon_form <- bf(PREDATION | trials(1) ~ SYMBIONT + (1 | BLOCK),
family = binomial(link = "logit")
)
# OR
mckeon_form <- bf(PREDATION ~ SYMBIONT + (1 | BLOCK),
family = bernoulli(link = "logit")
)
options(width = 150)
mckeon_form |> get_prior(data = mckeon) prior class coef group resp dpar nlpar lb ub source
(flat) b default
(flat) b SYMBIONTboth (vectorized)
(flat) b SYMBIONTcrabs (vectorized)
(flat) b SYMBIONTshrimp (vectorized)
student_t(3, 0, 2.5) Intercept default
student_t(3, 0, 2.5) sd 0 default
student_t(3, 0, 2.5) sd BLOCK 0 (vectorized)
student_t(3, 0, 2.5) sd Intercept BLOCK 0 (vectorized)options(width = 80)
mckeon |>
group_by(SYMBIONT) |>
summarise(
Mean = logit(mean(PREDATION)),
Median = logit(median(PREDATION)),
sd = logit(abs(sd(PREDATION)))
)# A tibble: 4 × 4
SYMBIONT Mean Median sd
<fct> <dbl> <dbl> <dbl>
1 none 2.20 Inf -0.810
2 crabs 0.405 Inf 0.0105
3 shrimp 0.201 Inf 0.0417
4 both -0.201 -Inf 0.0417 (Intercept) SYMBIONTcrabs SYMBIONTshrimp SYMBIONTboth
Inf 5.737305 5.737305 5.737305 mckeon |>
group_by(SYMBIONT) |>
summarise(
mean_response = mean(PREDATION),
mad_response = mad(PREDATION),
sd_response = sd(PREDATION),
) |>
mutate(
mean_logit = logit(mean_response),
# Delta method approximation
sd_logit = sd_response / (mean_response * (1 - mean_response))
)# A tibble: 4 × 6
SYMBIONT mean_response mad_response sd_response mean_logit sd_logit
<fct> <dbl> <dbl> <dbl> <dbl> <dbl>
1 none 0.9 0 0.308 2.20 3.42
2 crabs 0.6 0 0.503 0.405 2.09
3 shrimp 0.55 0 0.510 0.201 2.06
4 both 0.45 0 0.510 -0.201 2.06
mckeon_form <- bf(PREDATION | trials(1) ~ SYMBIONT + (1 | BLOCK),
family = binomial(link = "logit")
)
# OR
mckeon_form <- bf(PREDATION ~ SYMBIONT + (1 | BLOCK),
family = bernoulli(link = "logit")
)
get_prior(mckeon_form, data = mckeon) prior class coef group resp dpar nlpar lb ub
(flat) b
(flat) b SYMBIONTboth
(flat) b SYMBIONTcrabs
(flat) b SYMBIONTshrimp
student_t(3, 0, 2.5) Intercept
student_t(3, 0, 2.5) sd 0
student_t(3, 0, 2.5) sd BLOCK 0
student_t(3, 0, 2.5) sd Intercept BLOCK 0
source
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)priors <-
prior(normal(0, 2.5), class = "Intercept") +
prior(normal(0, 3), class = "b") +
prior(student_t(3, 0, 1.5), class = "sd")
priors <-
prior(normal(0, 2), class = "Intercept") +
prior(normal(0, 3), class = "b") +
prior(student_t(3, 0, 1.5), class = "sd")
mckeon_brm2 <- brm(mckeon_form,
data = mckeon,
prior = priors,
sample_prior = "only",
iter = 5000,
warmup = 1000,
chains = 3, cores = 3,
thin = 5,
control = list(adapt_delta = 0.99, max_treedepth = 20),
refresh = 0,
backend = "rstan"
)Compiling Stan program...Start sampling
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 failData points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
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 failData points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
mckeon_brm2 |>
emmeans(~SYMBIONT, type = "link") |>
as.data.frame() |>
ggplot(aes(y = emmean, x = SYMBIONT)) +
geom_hline(yintercept = c(5, -5), linetype = "dashed") +
geom_pointrange(aes(ymin = lower.HPD, ymax = upper.HPD)) +
geom_point(
data = mckeon, aes(y = PREDATION),
position = position_jitter(width = 0.2, height = 0),
alpha = 0.4, color = "red"
)Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
Now lets refit, conditioning on the data.
The desired updates require recompiling the modelCompiling Stan program...Start sampling
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 failData points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
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 failData points may overlap. Use the `jitter` argument to add some amount of
random variation to the location of data points and avoid overplotting.
mckeon_brm3 |>
emmeans(~SYMBIONT, type = "link") |>
as.data.frame() |>
ggplot(aes(y = emmean, x = SYMBIONT)) +
geom_hline(yintercept = c(5, -5), linetype = "dashed") +
geom_pointrange(aes(ymin = lower.HPD, ymax = upper.HPD)) +
geom_point(
data = mckeon, aes(y = PREDATION),
position = position_jitter(width = 0.2, height = 0),
alpha = 0.4, color = "red"
)Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.
[1] "b_Intercept" "b_SYMBIONTcrabs" "b_SYMBIONTshrimp"
[4] "b_SYMBIONTboth" "sd_BLOCK__Intercept" "Intercept"
[7] "r_BLOCK[1,Intercept]" "r_BLOCK[2,Intercept]" "r_BLOCK[3,Intercept]"
[10] "r_BLOCK[4,Intercept]" "r_BLOCK[5,Intercept]" "r_BLOCK[6,Intercept]"
[13] "r_BLOCK[7,Intercept]" "r_BLOCK[8,Intercept]" "r_BLOCK[9,Intercept]"
[16] "r_BLOCK[10,Intercept]" "prior_Intercept" "prior_b"
[19] "prior_sd_BLOCK" "lprior" "lp__"
[22] "accept_stat__" "stepsize__" "treedepth__"
[25] "n_leapfrog__" "divergent__" "energy__" 
While we are here, we might like to explore a random intercept/slope model
mckeon_form <- bf(PREDATION | trials(1) ~ SYMBIONT + (SYMBIONT | BLOCK),
family = binomial(link = "logit")
)
# OR
mckeon_form <- bf(PREDATION ~ SYMBIONT + (SYMBIONT | BLOCK),
family = bernoulli(link = "logit")
)
get_prior(mckeon_form, mckeon) prior class coef group resp dpar nlpar lb ub
(flat) b
(flat) b SYMBIONTboth
(flat) b SYMBIONTcrabs
(flat) b SYMBIONTshrimp
lkj(1) cor
lkj(1) cor BLOCK
student_t(3, 0, 2.5) Intercept
student_t(3, 0, 2.5) sd 0
student_t(3, 0, 2.5) sd BLOCK 0
student_t(3, 0, 2.5) sd Intercept BLOCK 0
student_t(3, 0, 2.5) sd SYMBIONTboth BLOCK 0
student_t(3, 0, 2.5) sd SYMBIONTcrabs BLOCK 0
student_t(3, 0, 2.5) sd SYMBIONTshrimp BLOCK 0
source
default
(vectorized)
(vectorized)
(vectorized)
default
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)## As there are not many observations, the following might be too ambitious without
## stronger priors
priors <-
prior(normal(0, 2), class = "Intercept") +
prior(normal(0, 3), class = "b") +
## prior(student_t(3,0,1.5), class = 'sd', coef = "Intercept", group = "BLOCK") +
prior(student_t(3, 0, 1.5), class = "sd") +
prior(lkj_corr_cholesky(1), class = "cor")
mckeon_brm4 <- brm(mckeon_form,
data = mckeon,
prior = priors,
sample_prior = "yes",
iter = 5000,
warmup = 1000,
chains = 3, cores = 3,
thin = 5,
refresh = 0,
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "rstan"
)Compiling Stan program...Start sampling [1] "b_Intercept"
[2] "b_SYMBIONTcrabs"
[3] "b_SYMBIONTshrimp"
[4] "b_SYMBIONTboth"
[5] "sd_BLOCK__Intercept"
[6] "sd_BLOCK__SYMBIONTcrabs"
[7] "sd_BLOCK__SYMBIONTshrimp"
[8] "sd_BLOCK__SYMBIONTboth"
[9] "cor_BLOCK__Intercept__SYMBIONTcrabs"
[10] "cor_BLOCK__Intercept__SYMBIONTshrimp"
[11] "cor_BLOCK__SYMBIONTcrabs__SYMBIONTshrimp"
[12] "cor_BLOCK__Intercept__SYMBIONTboth"
[13] "cor_BLOCK__SYMBIONTcrabs__SYMBIONTboth"
[14] "cor_BLOCK__SYMBIONTshrimp__SYMBIONTboth"
[15] "Intercept"
[16] "r_BLOCK[1,Intercept]"
[17] "r_BLOCK[2,Intercept]"
[18] "r_BLOCK[3,Intercept]"
[19] "r_BLOCK[4,Intercept]"
[20] "r_BLOCK[5,Intercept]"
[21] "r_BLOCK[6,Intercept]"
[22] "r_BLOCK[7,Intercept]"
[23] "r_BLOCK[8,Intercept]"
[24] "r_BLOCK[9,Intercept]"
[25] "r_BLOCK[10,Intercept]"
[26] "r_BLOCK[1,SYMBIONTcrabs]"
[27] "r_BLOCK[2,SYMBIONTcrabs]"
[28] "r_BLOCK[3,SYMBIONTcrabs]"
[29] "r_BLOCK[4,SYMBIONTcrabs]"
[30] "r_BLOCK[5,SYMBIONTcrabs]"
[31] "r_BLOCK[6,SYMBIONTcrabs]"
[32] "r_BLOCK[7,SYMBIONTcrabs]"
[33] "r_BLOCK[8,SYMBIONTcrabs]"
[34] "r_BLOCK[9,SYMBIONTcrabs]"
[35] "r_BLOCK[10,SYMBIONTcrabs]"
[36] "r_BLOCK[1,SYMBIONTshrimp]"
[37] "r_BLOCK[2,SYMBIONTshrimp]"
[38] "r_BLOCK[3,SYMBIONTshrimp]"
[39] "r_BLOCK[4,SYMBIONTshrimp]"
[40] "r_BLOCK[5,SYMBIONTshrimp]"
[41] "r_BLOCK[6,SYMBIONTshrimp]"
[42] "r_BLOCK[7,SYMBIONTshrimp]"
[43] "r_BLOCK[8,SYMBIONTshrimp]"
[44] "r_BLOCK[9,SYMBIONTshrimp]"
[45] "r_BLOCK[10,SYMBIONTshrimp]"
[46] "r_BLOCK[1,SYMBIONTboth]"
[47] "r_BLOCK[2,SYMBIONTboth]"
[48] "r_BLOCK[3,SYMBIONTboth]"
[49] "r_BLOCK[4,SYMBIONTboth]"
[50] "r_BLOCK[5,SYMBIONTboth]"
[51] "r_BLOCK[6,SYMBIONTboth]"
[52] "r_BLOCK[7,SYMBIONTboth]"
[53] "r_BLOCK[8,SYMBIONTboth]"
[54] "r_BLOCK[9,SYMBIONTboth]"
[55] "r_BLOCK[10,SYMBIONTboth]"
[56] "prior_Intercept"
[57] "prior_b"
[58] "prior_sd_BLOCK"
[59] "prior_cor_BLOCK"
[60] "lprior"
[61] "lp__"
[62] "accept_stat__"
[63] "stepsize__"
[64] "treedepth__"
[65] "n_leapfrog__"
[66] "divergent__"
[67] "energy__" 
Family: bernoulli
Links: mu = logit
Formula: PREDATION ~ SYMBIONT + (SYMBIONT | BLOCK)
Data: mckeon (Number of observations: 80)
Draws: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
total post-warmup draws = 2400
Multilevel Hyperparameters:
~BLOCK (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept) 2.15 1.19 0.20 4.87 1.00
sd(SYMBIONTcrabs) 3.40 3.47 0.12 11.58 1.00
sd(SYMBIONTshrimp) 2.39 2.24 0.09 8.24 1.00
sd(SYMBIONTboth) 2.41 2.33 0.08 8.11 1.00
cor(Intercept,SYMBIONTcrabs) 0.31 0.40 -0.55 0.91 1.00
cor(Intercept,SYMBIONTshrimp) 0.27 0.41 -0.66 0.88 1.00
cor(SYMBIONTcrabs,SYMBIONTshrimp) 0.38 0.41 -0.60 0.93 1.00
cor(Intercept,SYMBIONTboth) 0.21 0.41 -0.68 0.88 1.00
cor(SYMBIONTcrabs,SYMBIONTboth) 0.29 0.42 -0.68 0.90 1.00
cor(SYMBIONTshrimp,SYMBIONTboth) 0.34 0.42 -0.60 0.92 1.00
Bulk_ESS Tail_ESS
sd(Intercept) 1476 1361
sd(SYMBIONTcrabs) 1798 2290
sd(SYMBIONTshrimp) 1756 2160
sd(SYMBIONTboth) 1905 2208
cor(Intercept,SYMBIONTcrabs) 2144 2222
cor(Intercept,SYMBIONTshrimp) 2146 2089
cor(SYMBIONTcrabs,SYMBIONTshrimp) 2025 2110
cor(Intercept,SYMBIONTboth) 2198 2063
cor(SYMBIONTcrabs,SYMBIONTboth) 1804 2181
cor(SYMBIONTshrimp,SYMBIONTboth) 1606 2235
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.95 1.10 1.06 5.33 1.00 2169 2295
SYMBIONTcrabs -1.44 1.71 -4.31 2.42 1.00 2167 2133
SYMBIONTshrimp -2.19 1.49 -4.91 1.11 1.00 1843 2055
SYMBIONTboth -3.32 1.45 -6.32 -0.45 1.00 1795 2194
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).mckeon_brm4 |>
posterior_samples() |>
dplyr::select(-`lp__`) |>
pivot_longer(everything(), names_to = "key") |>
filter(!str_detect(key, "^r")) |>
mutate(
Type = ifelse(str_detect(key, "prior"), "Prior", "Posterior"),
## Class = ifelse(str_detect(key, 'Intercept'), 'Intercept',
## ifelse(str_detect(key, 'b'), 'b', 'sigma')),
Class = case_when(
str_detect(key, "(^b|^prior).*Intercept$") ~ "Intercept",
str_detect(key, "b_SYMBIONT.*|prior_b_SYMBIONT.*") &
!str_detect(key, ".*:.*") ~ "SYMBIONT",
str_detect(key, "sd") ~ "sd",
str_detect(key, "^cor|prior_cor") ~ "cor",
str_detect(key, "sigma") ~ "sigma"
),
Par = str_replace(key, "b_", "")
) |>
ggplot(aes(x = Type, y = value, color = Par)) +
stat_pointinterval(position = position_dodge()) +
facet_wrap(~Class, scales = "free")Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.
We can compare the two models using LOO
Computed from 2400 by 80 log-likelihood matrix.
Estimate SE
elpd_loo -27.8 6.9
p_loo 8.2 2.5
looic 55.6 13.8
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.Warning: Found 3 observations with a pareto_k > 0.7 in model 'mckeon_brm4'. We
recommend to set 'moment_match = TRUE' in order to perform moment matching for
problematic observations.
Computed from 2400 by 80 log-likelihood matrix.
Estimate SE
elpd_loo -25.9 6.1
p_loo 10.3 3.0
looic 51.7 12.2
------
MCSE of elpd_loo is NA.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
Pareto k diagnostic values:
Count Pct. Min. ESS
(-Inf, 0.7] (good) 77 96.2% 399
(0.7, 1] (bad) 3 3.8% <NA>
(1, Inf) (very bad) 0 0.0% <NA>
See help('pareto-k-diagnostic') for details. elpd_diff se_diff
mckeon_brm4 0.0 0.0
mckeon_brm3 -1.9 1.5 It appears that the random intercept/slope model is marginally better than the simpler random intercept model.
mckeon_brm4 %>%
posterior_samples() %>%
dplyr::select(-`lp__`) %>%
pivot_longer(everything(), names_to = "key") %>%
filter(!str_detect(key, "^r")) %>%
mutate(
Type = ifelse(str_detect(key, "prior"), "Prior", "Posterior"),
Class = case_when(
str_detect(key, "(^b|^prior).*Intercept$") ~ "Intercept",
str_detect(key, "b_SYMBIONT.*|prior_b") ~ "TREATMENT",
str_detect(key, "sd") ~ "sd",
str_detect(key, "^cor|prior_cor") ~ "cor",
str_detect(key, "sigma") ~ "sigma"
),
Par = str_replace(key, "b_", "")
) %>%
ggplot(aes(x = Type, y = value, color = Par)) +
stat_pointinterval(position = position_dodge(), show.legend = FALSE) +
facet_wrap(~Class, scales = "free")Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.
6 MCMC sampling diagnostics
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.
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_violinOf 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 <- mckeon_brm4 |> get_variables()
pars <- pars |>
str_extract("^b.Intercept|^b_SYMBIONT.*|[sS]igma|^sd.*") |>
na.omit()
pars[1] "b_Intercept" "b_SYMBIONTcrabs"
[3] "b_SYMBIONTshrimp" "b_SYMBIONTboth"
[5] "sd_BLOCK__Intercept" "sd_BLOCK__SYMBIONTcrabs"
[7] "sd_BLOCK__SYMBIONTshrimp" "sd_BLOCK__SYMBIONTboth"
attr(,"na.action")
[1] 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
[26] 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
[51] 59 60 61 62 63 64 65 66 67
attr(,"class")
[1] "omit"Warning: The following arguments were unrecognized and ignored: variablesNo divergences to plot.
# OR
mckeon_brm4 |> mcmc_plot(
type = "trace",
variable = "^b.Intercept|^b_SYMBIONT.*|[sS]igma|^sd.*",
regex = TRUE
)No divergences to plot.
The chains appear well mixed and very similar
- acf_bar (auto-correlation function): plots the auto-correlation between successive MCMC sample lags for each parameter and each chain
## OR
mckeon_brm4 |> mcmc_plot(
type = "acf_bar",
variable = "^b.Intercept|^b_SYMBIONT.*|[sS]igma|^sd.*",
regex = TRUE
)
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.
`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.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Ratios all very high.
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).
[1] "b_Intercept"
[2] "b_SYMBIONTcrabs"
[3] "b_SYMBIONTshrimp"
[4] "b_SYMBIONTboth"
[5] "sd_BLOCK__Intercept"
[6] "sd_BLOCK__SYMBIONTcrabs"
[7] "sd_BLOCK__SYMBIONTshrimp"
[8] "sd_BLOCK__SYMBIONTboth"
[9] "cor_BLOCK__Intercept__SYMBIONTcrabs"
[10] "cor_BLOCK__Intercept__SYMBIONTshrimp"
[11] "cor_BLOCK__SYMBIONTcrabs__SYMBIONTshrimp"
[12] "cor_BLOCK__Intercept__SYMBIONTboth"
[13] "cor_BLOCK__SYMBIONTcrabs__SYMBIONTboth"
[14] "cor_BLOCK__SYMBIONTshrimp__SYMBIONTboth"
[15] "Intercept"
[16] "r_BLOCK[1,Intercept]"
[17] "r_BLOCK[2,Intercept]"
[18] "r_BLOCK[3,Intercept]"
[19] "r_BLOCK[4,Intercept]"
[20] "r_BLOCK[5,Intercept]"
[21] "r_BLOCK[6,Intercept]"
[22] "r_BLOCK[7,Intercept]"
[23] "r_BLOCK[8,Intercept]"
[24] "r_BLOCK[9,Intercept]"
[25] "r_BLOCK[10,Intercept]"
[26] "r_BLOCK[1,SYMBIONTcrabs]"
[27] "r_BLOCK[2,SYMBIONTcrabs]"
[28] "r_BLOCK[3,SYMBIONTcrabs]"
[29] "r_BLOCK[4,SYMBIONTcrabs]"
[30] "r_BLOCK[5,SYMBIONTcrabs]"
[31] "r_BLOCK[6,SYMBIONTcrabs]"
[32] "r_BLOCK[7,SYMBIONTcrabs]"
[33] "r_BLOCK[8,SYMBIONTcrabs]"
[34] "r_BLOCK[9,SYMBIONTcrabs]"
[35] "r_BLOCK[10,SYMBIONTcrabs]"
[36] "r_BLOCK[1,SYMBIONTshrimp]"
[37] "r_BLOCK[2,SYMBIONTshrimp]"
[38] "r_BLOCK[3,SYMBIONTshrimp]"
[39] "r_BLOCK[4,SYMBIONTshrimp]"
[40] "r_BLOCK[5,SYMBIONTshrimp]"
[41] "r_BLOCK[6,SYMBIONTshrimp]"
[42] "r_BLOCK[7,SYMBIONTshrimp]"
[43] "r_BLOCK[8,SYMBIONTshrimp]"
[44] "r_BLOCK[9,SYMBIONTshrimp]"
[45] "r_BLOCK[10,SYMBIONTshrimp]"
[46] "r_BLOCK[1,SYMBIONTboth]"
[47] "r_BLOCK[2,SYMBIONTboth]"
[48] "r_BLOCK[3,SYMBIONTboth]"
[49] "r_BLOCK[4,SYMBIONTboth]"
[50] "r_BLOCK[5,SYMBIONTboth]"
[51] "r_BLOCK[6,SYMBIONTboth]"
[52] "r_BLOCK[7,SYMBIONTboth]"
[53] "r_BLOCK[8,SYMBIONTboth]"
[54] "r_BLOCK[9,SYMBIONTboth]"
[55] "r_BLOCK[10,SYMBIONTboth]"
[56] "prior_Intercept"
[57] "prior_b"
[58] "prior_sd_BLOCK"
[59] "prior_cor_BLOCK"
[60] "lprior"
[61] "lp__"
[62] "accept_stat__"
[63] "stepsize__"
[64] "treedepth__"
[65] "n_leapfrog__"
[66] "divergent__"
[67] "energy__" pars <- mckeon_brm4 |> get_variables()
pars <- str_extract(pars, "^b_.*|^sigma$|^sd.*") |> na.omit()
mckeon_brm4$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
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.
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.
Ratios all very high.
The ggmean package also as a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.
Of these, we will focus on:
- ggs_traceplot: 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).
The chains appear well mixed and very similar
- gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain
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.
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.
Ratios all very high.
7 Model validation
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.
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.
Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
instead.
The model draws appear to be consistent with the observed data.
- error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. Note, this is not really that useful for models that involve a binomial response
This is not really interpretable
- intervals: plots the observed data overlayed on top of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.
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.
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
mckeon_resids <- make_brms_dharma_res(mckeon_brm4, integerResponse = FALSE)
wrap_elements(~ testUniformity(mckeon_resids)) +
wrap_elements(~ plotResiduals(mckeon_resids, form = factor(rep(1, nrow(mckeon))))) +
wrap_elements(~ plotResiduals(mckeon_resids, quantreg = FALSE)) +
wrap_elements(~ testDispersion(mckeon_resids))
Conclusions:
- the simulated residuals do not suggest any issues with the residuals
- there is no evidence of a lack of fit.
8 Partial effects plots
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
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
## Partial residuals in binomial models are too confusing for the average viewer
## as they will yeild values that are not exactly 0 or 1 and this seems wrong.
## Partial.obs <- mckeon_brm3$data %>%
## mutate(Pred = predict(mckeon_brm3, re.form=NA)[,'Estimate'],
## Resid = resid(mckeon_brm3)[,'Estimate'],
## Obs = Pred + Resid)
mckeon_brm4 |>
epred_draws(newdata = mckeon, re_formula = NA) |>
median_hdci() |>
ggplot(aes(x = SYMBIONT, y = .epred)) +
geom_pointrange(aes(ymin = .lower, ymax = .upper)) +
geom_line() +
geom_point(
data = mckeon, aes(y = PREDATION, x = SYMBIONT),
alpha = 0.2,
position = position_jitter(width = 0.2, height = 0)
)
9 Model investigation
The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).
Family: bernoulli
Links: mu = logit
Formula: PREDATION ~ SYMBIONT + (SYMBIONT | BLOCK)
Data: mckeon (Number of observations: 80)
Draws: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
total post-warmup draws = 2400
Multilevel Hyperparameters:
~BLOCK (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept) 2.15 1.19 0.20 4.87 1.00
sd(SYMBIONTcrabs) 3.40 3.47 0.12 11.58 1.00
sd(SYMBIONTshrimp) 2.39 2.24 0.09 8.24 1.00
sd(SYMBIONTboth) 2.41 2.33 0.08 8.11 1.00
cor(Intercept,SYMBIONTcrabs) 0.31 0.40 -0.55 0.91 1.00
cor(Intercept,SYMBIONTshrimp) 0.27 0.41 -0.66 0.88 1.00
cor(SYMBIONTcrabs,SYMBIONTshrimp) 0.38 0.41 -0.60 0.93 1.00
cor(Intercept,SYMBIONTboth) 0.21 0.41 -0.68 0.88 1.00
cor(SYMBIONTcrabs,SYMBIONTboth) 0.29 0.42 -0.68 0.90 1.00
cor(SYMBIONTshrimp,SYMBIONTboth) 0.34 0.42 -0.60 0.92 1.00
Bulk_ESS Tail_ESS
sd(Intercept) 1476 1361
sd(SYMBIONTcrabs) 1798 2290
sd(SYMBIONTshrimp) 1756 2160
sd(SYMBIONTboth) 1905 2208
cor(Intercept,SYMBIONTcrabs) 2144 2222
cor(Intercept,SYMBIONTshrimp) 2146 2089
cor(SYMBIONTcrabs,SYMBIONTshrimp) 2025 2110
cor(Intercept,SYMBIONTboth) 2198 2063
cor(SYMBIONTcrabs,SYMBIONTboth) 1804 2181
cor(SYMBIONTshrimp,SYMBIONTboth) 1606 2235
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 2.95 1.10 1.06 5.33 1.00 2169 2295
SYMBIONTcrabs -1.44 1.71 -4.31 2.42 1.00 2167 2133
SYMBIONTshrimp -2.19 1.49 -4.91 1.11 1.00 1843 2055
SYMBIONTboth -3.32 1.45 -6.32 -0.45 1.00 1795 2194
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).Conclusions:
- the coefficients are presented on a logit scale. Whilst this is not relevant for the purpose of inference testing, it does make it difficult to interpret the coefficients.
- if we exponentiate the coefficients (\(log(\frac{\rho}{1-\rho})\) -> \(\frac{\rho}{1-\rho}\)), they will be presented on a odds ratio scale, and thus:
- the intercept (none symbionts) will be 19.12. That is, corals without a symbiont are 19.12 times more likely to be predated on than not predated on. The odds of predation in this the absence of symbionts is 19.12:1.
- in the presence of a crab symbiont, the odds of being predated on are only 0.24 times that of the none symbiont group. That is, in the presence of a crab symbiont, the odds of predation decline by 76%.
- in the presence of a shrimp symbiont, the odds of being predated on are only 0.11 times that of the none symbiont group. That is, in the presence of a shrimp symbiont, the odds of predation decline by 89%.
- in the presence of both crab and shrimp symbionts, the odds of being predated on are only 0.04 times that of the none symbiont group. That is, in the presence of both crab and shrimp symbiont, the odds of predation decline by 96%.
- if we back-transform the intercept full to the response scale (probability scale), (\(log(\frac{\rho}{1-\rho})\) -> \(\rho\)), the intercept is interpreted as the probability that corals will be predated in the absence of of symbionts is 1
# A draws_df: 800 iterations, 3 chains, and 61 variables
b_Intercept b_SYMBIONTcrabs b_SYMBIONTshrimp b_SYMBIONTboth
1 2.6 0.56 -1.52 -4.2
2 3.0 -3.58 -2.42 -2.8
3 2.1 -0.42 -0.37 -2.8
4 2.7 -1.13 -2.36 -1.7
5 2.1 1.07 -2.01 -1.7
6 3.3 -2.98 -2.56 -4.2
7 2.4 -1.62 -0.90 -1.9
8 2.5 -1.37 -2.22 -5.1
9 1.8 -0.38 -0.83 -2.3
10 2.6 -2.85 -3.29 -2.0
sd_BLOCK__Intercept sd_BLOCK__SYMBIONTcrabs sd_BLOCK__SYMBIONTshrimp
1 1.57 2.01 5.504
2 2.05 1.81 1.261
3 1.78 1.16 5.007
4 1.15 0.68 5.866
5 1.27 1.75 4.340
6 1.73 0.51 0.534
7 1.65 2.19 0.792
8 0.84 1.66 1.465
9 2.65 1.37 0.027
10 2.75 1.66 0.685
sd_BLOCK__SYMBIONTboth
1 1.062
2 3.219
3 3.351
4 5.368
5 4.295
6 0.725
7 0.049
8 4.671
9 1.176
10 0.349
# ... with 2390 more draws, and 53 more variables
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}| variable | median | lower | upper | Pl | Pg | rhat | ess_bulk | ess_tail |
|---|---|---|---|---|---|---|---|---|
| b_Intercept | 16.9795558 | 0.9493822 | 1.356772e+02 | 0.0016667 | 0.9983333 | 1.0005945 | 2169.079 | 2294.643 |
| b_SYMBIONTcrabs | 0.2005002 | 0.0012488 | 5.768128e+00 | 0.8150000 | 0.1850000 | 1.0007817 | 2167.227 | 2133.030 |
| b_SYMBIONTshrimp | 0.1009254 | 0.0009651 | 1.495900e+00 | 0.9233333 | 0.0766667 | 1.0016047 | 1842.846 | 2055.451 |
| b_SYMBIONTboth | 0.0368083 | 0.0001099 | 3.676436e-01 | 0.9858333 | 0.0141667 | 0.9996782 | 1796.009 | 2194.054 |
| sd_BLOCK__Intercept | 7.8170801 | 1.0024454 | 6.727319e+01 | 0.0000000 | 1.0000000 | 0.9998929 | 1475.617 | 1361.360 |
| sd_BLOCK__SYMBIONTcrabs | 12.9512616 | 1.0022307 | 8.953066e+03 | 0.0000000 | 1.0000000 | 0.9997778 | 1797.154 | 2290.099 |
| sd_BLOCK__SYMBIONTshrimp | 6.2784206 | 1.0028237 | 7.967616e+02 | 0.0000000 | 1.0000000 | 1.0017262 | 1756.293 | 2160.197 |
| sd_BLOCK__SYMBIONTboth | 6.2137576 | 1.0001731 | 7.487175e+02 | 0.0000000 | 1.0000000 | 1.0003384 | 1905.228 | 2207.855 |
| cor_BLOCK__Intercept__SYMBIONTcrabs | 1.4314022 | 0.5641743 | 2.463366e+00 | 0.2216667 | 0.7783333 | 1.0005094 | 2144.390 | 2222.299 |
| cor_BLOCK__Intercept__SYMBIONTshrimp | 1.3875938 | 0.4760768 | 2.356132e+00 | 0.2504167 | 0.7495833 | 0.9999908 | 2146.956 | 2088.625 |
| cor_BLOCK__SYMBIONTcrabs__SYMBIONTshrimp | 1.6052128 | 0.5468642 | 2.524864e+00 | 0.1858333 | 0.8141667 | 1.0011434 | 2025.956 | 2110.236 |
| cor_BLOCK__Intercept__SYMBIONTboth | 1.2788376 | 0.4328817 | 2.274023e+00 | 0.3083333 | 0.6916667 | 1.0001088 | 2198.585 | 2062.885 |
| cor_BLOCK__SYMBIONTcrabs__SYMBIONTboth | 1.4287933 | 0.4537353 | 2.367508e+00 | 0.2466667 | 0.7533333 | 0.9999072 | 1803.095 | 2181.215 |
| cor_BLOCK__SYMBIONTshrimp__SYMBIONTboth | 1.5295511 | 0.5726949 | 2.539087e+00 | 0.2133333 | 0.7866667 | 1.0005362 | 1605.541 | 2234.621 |
| Intercept | 3.3329118 | 0.2338775 | 1.731357e+01 | 0.0983333 | 0.9016667 | 1.0015745 | 2346.332 | 2187.828 |
| r_BLOCK[1,Intercept] | 0.0931594 | 0.0001148 | 9.200035e-01 | 0.9641667 | 0.0358333 | 1.0004264 | 1608.804 | 1478.134 |
| r_BLOCK[2,Intercept] | 0.2213524 | 0.0003908 | 1.610838e+00 | 0.8779167 | 0.1220833 | 1.0005555 | 1724.775 | 1976.003 |
| r_BLOCK[3,Intercept] | 0.2298054 | 0.0007686 | 1.469713e+00 | 0.8808333 | 0.1191667 | 0.9997249 | 1904.461 | 2120.734 |
| r_BLOCK[4,Intercept] | 0.2289850 | 0.0006309 | 1.556997e+00 | 0.8833333 | 0.1166667 | 1.0002289 | 1704.193 | 2289.127 |
| r_BLOCK[5,Intercept] | 0.9396625 | 0.0090167 | 6.843036e+00 | 0.5325000 | 0.4675000 | 0.9994026 | 2027.899 | 2411.620 |
| r_BLOCK[6,Intercept] | 2.5332838 | 0.0470625 | 4.466575e+01 | 0.1958333 | 0.8041667 | 0.9996394 | 2031.548 | 2250.614 |
| r_BLOCK[7,Intercept] | 7.1884377 | 0.0915078 | 4.584838e+02 | 0.0804167 | 0.9195833 | 1.0016462 | 1745.953 | 1871.426 |
| r_BLOCK[8,Intercept] | 7.1316871 | 0.0961562 | 5.819502e+02 | 0.0900000 | 0.9100000 | 1.0004001 | 1705.950 | 2191.186 |
| r_BLOCK[9,Intercept] | 6.9060295 | 0.0369566 | 7.324657e+02 | 0.0862500 | 0.9137500 | 1.0020705 | 1863.148 | 2197.916 |
| r_BLOCK[10,Intercept] | 1.8009988 | 0.0083170 | 2.378800e+01 | 0.2850000 | 0.7150000 | 1.0013829 | 1723.623 | 2095.223 |
| r_BLOCK[1,SYMBIONTcrabs] | 0.1057067 | 0.0000000 | 1.847332e+00 | 0.8595833 | 0.1404167 | 0.9999083 | 1696.087 | 2097.514 |
| r_BLOCK[2,SYMBIONTcrabs] | 0.1011589 | 0.0000000 | 1.621368e+00 | 0.8712500 | 0.1287500 | 1.0010480 | 1829.851 | 2137.464 |
| r_BLOCK[3,SYMBIONTcrabs] | 0.0947467 | 0.0000000 | 1.492297e+00 | 0.8758333 | 0.1241667 | 0.9996350 | 1581.459 | 1983.190 |
| r_BLOCK[4,SYMBIONTcrabs] | 0.0906858 | 0.0000000 | 1.505313e+00 | 0.8691667 | 0.1308333 | 1.0002267 | 1745.715 | 2377.100 |
| r_BLOCK[5,SYMBIONTcrabs] | 1.7407549 | 0.0050246 | 4.838996e+02 | 0.3258333 | 0.6741667 | 1.0030566 | 2259.308 | 1916.344 |
| r_BLOCK[6,SYMBIONTcrabs] | 2.6594594 | 0.0023820 | 5.718778e+03 | 0.2529167 | 0.7470833 | 0.9990437 | 2234.653 | 2113.964 |
| r_BLOCK[7,SYMBIONTcrabs] | 4.4115529 | 0.0001198 | 8.030981e+04 | 0.2258333 | 0.7741667 | 1.0003395 | 2060.188 | 2444.366 |
| r_BLOCK[8,SYMBIONTcrabs] | 4.3078418 | 0.0006814 | 6.842921e+04 | 0.2320833 | 0.7679167 | 1.0016648 | 2126.383 | 2369.811 |
| r_BLOCK[9,SYMBIONTcrabs] | 4.2914666 | 0.0017548 | 5.903314e+04 | 0.2216667 | 0.7783333 | 0.9997776 | 2064.066 | 2276.536 |
| r_BLOCK[10,SYMBIONTcrabs] | 4.0749951 | 0.0052980 | 4.549004e+04 | 0.2037500 | 0.7962500 | 1.0023311 | 1996.097 | 2189.067 |
| r_BLOCK[1,SYMBIONTshrimp] | 0.2391954 | 0.0000000 | 1.955260e+00 | 0.8158333 | 0.1841667 | 1.0009367 | 1914.814 | 2078.217 |
| r_BLOCK[2,SYMBIONTshrimp] | 0.2436122 | 0.0000000 | 1.696204e+00 | 0.8258333 | 0.1741667 | 1.0017388 | 1614.421 | 2310.313 |
| r_BLOCK[3,SYMBIONTshrimp] | 0.2170023 | 0.0000000 | 1.633490e+00 | 0.8470833 | 0.1529167 | 1.0019608 | 1813.510 | 2393.158 |
| r_BLOCK[4,SYMBIONTshrimp] | 0.2273558 | 0.0000000 | 1.733808e+00 | 0.8245833 | 0.1754167 | 1.0027234 | 1683.151 | 2146.537 |
| r_BLOCK[5,SYMBIONTshrimp] | 0.8320058 | 0.0000875 | 6.242993e+00 | 0.6029167 | 0.3970833 | 1.0006494 | 2311.938 | 2251.069 |
| r_BLOCK[6,SYMBIONTshrimp] | 1.9216520 | 0.0013234 | 3.617197e+02 | 0.2533333 | 0.7466667 | 1.0008797 | 2294.319 | 2257.701 |
| r_BLOCK[7,SYMBIONTshrimp] | 2.7243959 | 0.0009955 | 3.204963e+03 | 0.2404167 | 0.7595833 | 1.0009013 | 2186.625 | 2243.150 |
| r_BLOCK[8,SYMBIONTshrimp] | 2.6384816 | 0.0011415 | 6.302724e+03 | 0.2220833 | 0.7779167 | 1.0006606 | 2038.689 | 2448.646 |
| r_BLOCK[9,SYMBIONTshrimp] | 2.7451160 | 0.0081512 | 2.616071e+03 | 0.2329167 | 0.7670833 | 1.0014791 | 2011.952 | 2253.266 |
| r_BLOCK[10,SYMBIONTshrimp] | 3.0572310 | 0.0189161 | 2.080085e+03 | 0.1991667 | 0.8008333 | 1.0012863 | 1866.872 | 2266.207 |
| r_BLOCK[1,SYMBIONTboth] | 0.3379118 | 0.0000000 | 3.223799e+00 | 0.7625000 | 0.2375000 | 1.0007930 | 2051.365 | 2329.975 |
| r_BLOCK[2,SYMBIONTboth] | 0.3256919 | 0.0000000 | 2.546152e+00 | 0.7962500 | 0.2037500 | 1.0009283 | 1874.095 | 2117.824 |
| r_BLOCK[3,SYMBIONTboth] | 0.3242087 | 0.0000000 | 2.416801e+00 | 0.7937500 | 0.2062500 | 1.0003547 | 1742.901 | 1884.533 |
| r_BLOCK[4,SYMBIONTboth] | 0.3265888 | 0.0000000 | 2.329331e+00 | 0.7970833 | 0.2029167 | 1.0005674 | 1740.111 | 2288.133 |
| r_BLOCK[5,SYMBIONTboth] | 0.4593611 | 0.0000000 | 2.512963e+00 | 0.7495833 | 0.2504167 | 1.0007437 | 2093.274 | 2447.752 |
| r_BLOCK[6,SYMBIONTboth] | 1.0142797 | 0.0000025 | 1.312239e+01 | 0.4845833 | 0.5154167 | 0.9994514 | 2357.973 | 2187.028 |
| r_BLOCK[7,SYMBIONTboth] | 2.7863000 | 0.0001577 | 2.434989e+03 | 0.2362500 | 0.7637500 | 1.0006519 | 2150.135 | 2028.503 |
| r_BLOCK[8,SYMBIONTboth] | 2.9442757 | 0.0002624 | 2.369287e+03 | 0.2233333 | 0.7766667 | 1.0020189 | 1892.232 | 2367.680 |
| r_BLOCK[9,SYMBIONTboth] | 3.1410159 | 0.0008290 | 2.045463e+03 | 0.2120833 | 0.7879167 | 1.0020604 | 2252.120 | 2370.097 |
| r_BLOCK[10,SYMBIONTboth] | 3.9591838 | 0.0165197 | 1.895715e+03 | 0.1704167 | 0.8295833 | 1.0005154 | 1845.681 | 2188.774 |
| prior_Intercept | 0.9683435 | 0.0010802 | 2.408192e+01 | 0.5037500 | 0.4962500 | 0.9998291 | 2244.592 | 2189.815 |
| prior_b | 0.9420457 | 0.0000864 | 1.430789e+02 | 0.5083333 | 0.4916667 | 1.0010273 | 2531.765 | 2436.551 |
| prior_sd_BLOCK | 3.0820146 | 1.0001254 | 9.197562e+01 | 0.0000000 | 1.0000000 | 0.9998882 | 2310.486 | 1987.187 |
| prior_cor_BLOCK | 1.0038141 | 0.3930119 | 2.057522e+00 | 0.4962500 | 0.5037500 | 1.0004393 | 2298.553 | 1978.868 |
| lprior | 0.0000000 | 0.0000000 | 0.000000e+00 | 1.0000000 | 0.0000000 | 0.9999232 | 1803.622 | 2310.997 |
| lp__ | 0.0000000 | 0.0000000 | 0.000000e+00 | 1.0000000 | 0.0000000 | 0.9998792 | 1875.493 | 2267.273 |
Warning: Dropping 'draws_df' class as required metadata was removed.| variable | median | lower | upper | Pl | Pg | rhat | ess_bulk | ess_tail |
|---|---|---|---|---|---|---|---|---|
| b_Intercept | 16.9795558 | 0.9493822 | 135.6772240 | 0.0016667 | 0.9983333 | 1.0007811 | 2159.224 | 2254.273 |
| b_SYMBIONTcrabs | 0.2005002 | 0.0012488 | 5.7681283 | 0.8150000 | 0.1850000 | 1.0002414 | 2151.557 | 2107.985 |
| b_SYMBIONTshrimp | 0.1009254 | 0.0009651 | 1.4959003 | 0.9233333 | 0.0766667 | 1.0030755 | 1783.447 | 2024.067 |
| b_SYMBIONTboth | 0.0368083 | 0.0001099 | 0.3676436 | 0.9858333 | 0.0141667 | 1.0000046 | 1763.825 | 2166.855 |
| sd_BLOCK__Intercept | 7.8170801 | 1.0024454 | 67.2731901 | 0.0000000 | 1.0000000 | 1.0002783 | 1467.218 | 1354.369 |
| sd_BLOCK__SYMBIONTcrabs | 12.9512616 | 1.0022307 | 8953.0662531 | 0.0000000 | 1.0000000 | 0.9997799 | 1768.058 | 2255.879 |
| sd_BLOCK__SYMBIONTshrimp | 6.2784206 | 1.0028237 | 796.7615563 | 0.0000000 | 1.0000000 | 1.0026563 | 1721.165 | 2143.960 |
| sd_BLOCK__SYMBIONTboth | 6.2137576 | 1.0001731 | 748.7175253 | 0.0000000 | 1.0000000 | 1.0007922 | 1885.999 | 2195.443 |
mckeon_brm4$fit |>
tidyMCMC(
estimate.method = "median",
conf.int = TRUE, conf.method = "HPDinterval",
rhat = TRUE, ess = TRUE
)# A tibble: 60 × 7
term estimate std.error conf.low conf.high rhat ess
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
1 b_Intercept 2.95 1.10 8.82e-1 5.15 1.000 2173
2 b_SYMBIONTcrabs -1.44 1.71 -4.24e+0 2.46 1.000 2083
3 b_SYMBIONTshrimp -2.19 1.49 -5.33e+0 0.655 1.00 1822
4 b_SYMBIONTboth -3.32 1.45 -6.26e+0 -0.419 1.000 1815
5 sd_BLOCK__Intercept 2.15 1.19 2.44e-3 4.21 1.000 1634
6 sd_BLOCK__SYMBIONTcrabs 3.40 3.47 2.23e-3 9.10 0.999 1871
7 sd_BLOCK__SYMBIONTshrimp 2.39 2.24 2.82e-3 6.68 1.00 1959
8 sd_BLOCK__SYMBIONTboth 2.41 2.33 1.73e-4 6.62 0.999 2097
9 cor_BLOCK__Intercept__SYMB… 0.306 0.396 -4.52e-1 0.945 0.999 2149
10 cor_BLOCK__Intercept__SYMB… 0.271 0.409 -5.11e-1 0.934 1.000 2162
# ℹ 50 more rowsConclusions:
see above
Due to the presence of a log transform in the predictor, it is better to use the regex version.
[1] "b_Intercept"
[2] "b_SYMBIONTcrabs"
[3] "b_SYMBIONTshrimp"
[4] "b_SYMBIONTboth"
[5] "sd_BLOCK__Intercept"
[6] "sd_BLOCK__SYMBIONTcrabs"
[7] "sd_BLOCK__SYMBIONTshrimp"
[8] "sd_BLOCK__SYMBIONTboth"
[9] "cor_BLOCK__Intercept__SYMBIONTcrabs"
[10] "cor_BLOCK__Intercept__SYMBIONTshrimp"
[11] "cor_BLOCK__SYMBIONTcrabs__SYMBIONTshrimp"
[12] "cor_BLOCK__Intercept__SYMBIONTboth"
[13] "cor_BLOCK__SYMBIONTcrabs__SYMBIONTboth"
[14] "cor_BLOCK__SYMBIONTshrimp__SYMBIONTboth"
[15] "Intercept"
[16] "r_BLOCK[1,Intercept]"
[17] "r_BLOCK[2,Intercept]"
[18] "r_BLOCK[3,Intercept]"
[19] "r_BLOCK[4,Intercept]"
[20] "r_BLOCK[5,Intercept]"
[21] "r_BLOCK[6,Intercept]"
[22] "r_BLOCK[7,Intercept]"
[23] "r_BLOCK[8,Intercept]"
[24] "r_BLOCK[9,Intercept]"
[25] "r_BLOCK[10,Intercept]"
[26] "r_BLOCK[1,SYMBIONTcrabs]"
[27] "r_BLOCK[2,SYMBIONTcrabs]"
[28] "r_BLOCK[3,SYMBIONTcrabs]"
[29] "r_BLOCK[4,SYMBIONTcrabs]"
[30] "r_BLOCK[5,SYMBIONTcrabs]"
[31] "r_BLOCK[6,SYMBIONTcrabs]"
[32] "r_BLOCK[7,SYMBIONTcrabs]"
[33] "r_BLOCK[8,SYMBIONTcrabs]"
[34] "r_BLOCK[9,SYMBIONTcrabs]"
[35] "r_BLOCK[10,SYMBIONTcrabs]"
[36] "r_BLOCK[1,SYMBIONTshrimp]"
[37] "r_BLOCK[2,SYMBIONTshrimp]"
[38] "r_BLOCK[3,SYMBIONTshrimp]"
[39] "r_BLOCK[4,SYMBIONTshrimp]"
[40] "r_BLOCK[5,SYMBIONTshrimp]"
[41] "r_BLOCK[6,SYMBIONTshrimp]"
[42] "r_BLOCK[7,SYMBIONTshrimp]"
[43] "r_BLOCK[8,SYMBIONTshrimp]"
[44] "r_BLOCK[9,SYMBIONTshrimp]"
[45] "r_BLOCK[10,SYMBIONTshrimp]"
[46] "r_BLOCK[1,SYMBIONTboth]"
[47] "r_BLOCK[2,SYMBIONTboth]"
[48] "r_BLOCK[3,SYMBIONTboth]"
[49] "r_BLOCK[4,SYMBIONTboth]"
[50] "r_BLOCK[5,SYMBIONTboth]"
[51] "r_BLOCK[6,SYMBIONTboth]"
[52] "r_BLOCK[7,SYMBIONTboth]"
[53] "r_BLOCK[8,SYMBIONTboth]"
[54] "r_BLOCK[9,SYMBIONTboth]"
[55] "r_BLOCK[10,SYMBIONTboth]"
[56] "prior_Intercept"
[57] "prior_b"
[58] "prior_sd_BLOCK"
[59] "prior_cor_BLOCK"
[60] "lprior"
[61] "lp__"
[62] "accept_stat__"
[63] "stepsize__"
[64] "treedepth__"
[65] "n_leapfrog__"
[66] "divergent__"
[67] "energy__" # A tibble: 9,600 × 5
# Groups: .variable [4]
.chain .iteration .draw .variable .value
<int> <int> <int> <chr> <dbl>
1 1 1 1 b_Intercept 2.61
2 1 2 2 b_Intercept 2.98
3 1 3 3 b_Intercept 2.14
4 1 4 4 b_Intercept 2.68
5 1 5 5 b_Intercept 2.09
6 1 6 6 b_Intercept 3.29
7 1 7 7 b_Intercept 2.41
8 1 8 8 b_Intercept 2.50
9 1 9 9 b_Intercept 1.78
10 1 10 10 b_Intercept 2.60
# ℹ 9,590 more rowsWe can then summarise this
# A tibble: 4 × 7
.variable .value .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 b_Intercept 2.83 0.882 5.15 0.95 median hdci
2 b_SYMBIONTboth -3.30 -6.19 -0.335 0.95 median hdci
3 b_SYMBIONTcrabs -1.61 -4.24 2.46 0.95 median hdci
4 b_SYMBIONTshrimp -2.29 -5.15 0.857 0.95 median hdci # A tibble: 4 × 7
.variable .value .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 b_Intercept 17.0 0.949 136. 0.95 median hdci
2 b_SYMBIONTboth 0.0368 0.000110 0.367 0.95 median hdci
3 b_SYMBIONTcrabs 0.201 0.00125 5.74 0.95 median hdci
4 b_SYMBIONTshrimp 0.101 0.000965 1.48 0.95 median hdci mckeon_brm4 |>
gather_draws(`b_Intercept.*|b_SYMBIONT.*`, regex = TRUE) |>
ggplot() +
geom_vline(xintercept = 0, linetype = "dashed") +
stat_slab(aes(
x = .value, y = .variable,
fill = stat(ggdist::cut_cdf_qi(cdf,
.width = c(0.5, 0.8, 0.95),
labels = scales::percent_format()
))
), color = "black") +
scale_fill_brewer("Interval", direction = -1, na.translate = FALSE)Warning: `stat(ggdist::cut_cdf_qi(cdf, .width = c(0.5, 0.8, 0.95), labels =
scales::percent_format()))` was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(ggdist::cut_cdf_qi(cdf, .width = c(0.5, 0.8, 0.95),
labels = scales::percent_format()))` instead.
mckeon_brm4 |>
gather_draws(`.Intercept.*|b_SYMBIONT.*`, regex = TRUE) |>
ggplot() +
geom_vline(xintercept = 0, linetype = "dashed") +
stat_halfeye(aes(x = .value, y = .variable)) +
theme_classic()
mckeon_brm4 |>
gather_draws(`^b_.*`, regex = TRUE) |>
filter(.variable != "b_Intercept") |>
ggplot() +
stat_halfeye(aes(x = .value, y = .variable)) +
facet_wrap(~.variable, scales = "free")
mckeon_brm4 |>
gather_draws(`^b_.*`, regex = TRUE) |>
filter(.variable != "b_Intercept") |>
ggplot() +
stat_halfeye(aes(x = .value, y = .variable)) +
geom_vline(xintercept = 0, linetype = "dashed")
mckeon_brm4 |>
gather_draws(`^b_.*`, regex = TRUE) |>
filter(.variable != "b_Intercept") |>
ggplot() +
geom_density_ridges(aes(x = .value, y = .variable), alpha = 0.4) +
geom_vline(xintercept = 0, linetype = "dashed")Picking joint bandwidth of 0.278
## Or in colour
mckeon_brm4 |>
gather_draws(`^b_.*`, regex = TRUE) |>
filter(.variable != "b_Intercept") |>
ggplot() +
geom_density_ridges_gradient(
aes(
x = (.value),
y = .variable,
fill = stat(x)
),
alpha = 0.4, colour = "white",
quantile_lines = TRUE,
quantiles = c(0.025, 0.975)
) +
geom_vline(xintercept = 1, linetype = "dashed") +
scale_x_continuous() +
scale_fill_viridis_c(option = "C")Picking joint bandwidth of 0.278
## Fractional scale
mckeon_brm4 |>
gather_draws(`^b_.*`, regex = TRUE) |>
filter(.variable != "b_Intercept") |>
ggplot() +
geom_density_ridges_gradient(
aes(
x = exp(.value),
y = .variable,
fill = stat(x)
),
alpha = 0.4, colour = "white",
quantile_lines = TRUE,
quantiles = c(0.025, 0.975)
) +
geom_vline(xintercept = 1, linetype = "dashed") +
scale_x_continuous(trans = scales::log2_trans()) +
scale_fill_viridis_c(option = "C")Picking joint bandwidth of 0.401
# A tibble: 2,400 × 70
.chain .iteration .draw b_Intercept b_SYMBIONTcrabs b_SYMBIONTshrimp
<int> <int> <int> <dbl> <dbl> <dbl>
1 1 1 1 2.61 0.558 -1.52
2 1 2 2 2.98 -3.58 -2.42
3 1 3 3 2.14 -0.416 -0.366
4 1 4 4 2.68 -1.13 -2.36
5 1 5 5 2.09 1.07 -2.01
6 1 6 6 3.29 -2.98 -2.56
7 1 7 7 2.41 -1.62 -0.897
8 1 8 8 2.50 -1.37 -2.22
9 1 9 9 1.78 -0.378 -0.835
10 1 10 10 2.60 -2.85 -3.29
# ℹ 2,390 more rows
# ℹ 64 more variables: b_SYMBIONTboth <dbl>, sd_BLOCK__Intercept <dbl>,
# sd_BLOCK__SYMBIONTcrabs <dbl>, sd_BLOCK__SYMBIONTshrimp <dbl>,
# sd_BLOCK__SYMBIONTboth <dbl>, cor_BLOCK__Intercept__SYMBIONTcrabs <dbl>,
# cor_BLOCK__Intercept__SYMBIONTshrimp <dbl>,
# cor_BLOCK__SYMBIONTcrabs__SYMBIONTshrimp <dbl>,
# cor_BLOCK__Intercept__SYMBIONTboth <dbl>, …# A tibble: 2,400 × 23
.chain .iteration .draw b_Intercept b_SYMBIONTcrabs b_SYMBIONTshrimp
<int> <int> <int> <dbl> <dbl> <dbl>
1 1 1 1 2.61 0.558 -1.52
2 1 2 2 2.98 -3.58 -2.42
3 1 3 3 2.14 -0.416 -0.366
4 1 4 4 2.68 -1.13 -2.36
5 1 5 5 2.09 1.07 -2.01
6 1 6 6 3.29 -2.98 -2.56
7 1 7 7 2.41 -1.62 -0.897
8 1 8 8 2.50 -1.37 -2.22
9 1 9 9 1.78 -0.378 -0.835
10 1 10 10 2.60 -2.85 -3.29
# ℹ 2,390 more rows
# ℹ 17 more variables: b_SYMBIONTboth <dbl>, sd_BLOCK__Intercept <dbl>,
# cor_BLOCK__Intercept__SYMBIONTcrabs <dbl>,
# cor_BLOCK__Intercept__SYMBIONTshrimp <dbl>,
# cor_BLOCK__Intercept__SYMBIONTboth <dbl>, Intercept <dbl>,
# `r_BLOCK[1,Intercept]` <dbl>, `r_BLOCK[2,Intercept]` <dbl>,
# `r_BLOCK[3,Intercept]` <dbl>, `r_BLOCK[4,Intercept]` <dbl>, …Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.# A tibble: 2,400 × 61
b_Intercept b_SYMBIONTcrabs b_SYMBIONTshrimp b_SYMBIONTboth
<dbl> <dbl> <dbl> <dbl>
1 2.61 0.558 -1.52 -4.24
2 2.98 -3.58 -2.42 -2.78
3 2.14 -0.416 -0.366 -2.76
4 2.68 -1.13 -2.36 -1.66
5 2.09 1.07 -2.01 -1.70
6 3.29 -2.98 -2.56 -4.25
7 2.41 -1.62 -0.897 -1.92
8 2.50 -1.37 -2.22 -5.08
9 1.78 -0.378 -0.835 -2.31
10 2.60 -2.85 -3.29 -1.95
# ℹ 2,390 more rows
# ℹ 57 more variables: sd_BLOCK__Intercept <dbl>,
# sd_BLOCK__SYMBIONTcrabs <dbl>, sd_BLOCK__SYMBIONTshrimp <dbl>,
# sd_BLOCK__SYMBIONTboth <dbl>, cor_BLOCK__Intercept__SYMBIONTcrabs <dbl>,
# cor_BLOCK__Intercept__SYMBIONTshrimp <dbl>,
# cor_BLOCK__SYMBIONTcrabs__SYMBIONTshrimp <dbl>,
# cor_BLOCK__Intercept__SYMBIONTboth <dbl>, … y ymin ymax .width .point .interval
1 0.1798171 0.008074149 0.321401 0.95 median hdci y ymin ymax .width .point .interval
1 0.4706625 0.1406366 0.734526 0.95 median hdci## if we had random intercept/slope
mckeon_brm4 |>
bayes_R2(re.form = ~ (SYMBIONT | BLOCK), summary = FALSE) |>
median_hdci() y ymin ymax .width .point .interval
1 0.7080777 0.5762833 0.8229876 0.95 median hdcimckeon_brm4 |> modelsummary(
statistic = c("conf.low", "conf.high"),
shape = term ~ statistic,
exponentiate = FALSE
)Warning:
`modelsummary` uses the `performance` package to extract goodness-of-fit
statistics from models of this class. You can specify the statistics you wish
to compute by supplying a `metrics` argument to `modelsummary`, which will then
push it forward to `performance`. Acceptable values are: "all", "common",
"none", or a character vector of metrics names. For example: `modelsummary(mod,
metrics = c("RMSE", "R2")` Note that some metrics are computationally
expensive. See `?performance::performance` for details.
This warning appears once per session.| (1) | |||
|---|---|---|---|
| Est. | 2.5 % | 97.5 % | |
| b_Intercept | 2.832 | 1.056 | 5.327 |
| b_SYMBIONTcrabs | -1.607 | -4.309 | 2.421 |
| b_SYMBIONTshrimp | -2.293 | -4.908 | 1.110 |
| b_SYMBIONTboth | -3.302 | -6.322 | -0.447 |
| sd_BLOCK__Intercept | 2.056 | 0.204 | 4.865 |
| sd_BLOCK__SYMBIONTcrabs | 2.561 | 0.119 | 11.582 |
| sd_BLOCK__SYMBIONTshrimp | 1.837 | 0.088 | 8.240 |
| sd_BLOCK__SYMBIONTboth | 1.827 | 0.081 | 8.105 |
| cor_BLOCK__Intercept__SYMBIONTcrabs | 0.359 | -0.554 | 0.907 |
| cor_BLOCK__Intercept__SYMBIONTshrimp | 0.328 | -0.662 | 0.882 |
| cor_BLOCK__SYMBIONTcrabs__SYMBIONTshrimp | 0.473 | -0.602 | 0.928 |
| cor_BLOCK__Intercept__SYMBIONTboth | 0.246 | -0.682 | 0.877 |
| cor_BLOCK__SYMBIONTcrabs__SYMBIONTboth | 0.357 | -0.677 | 0.900 |
| cor_BLOCK__SYMBIONTshrimp__SYMBIONTboth | 0.425 | -0.604 | 0.921 |
| Num.Obs. | 80 | ||
| R2 | 0.708 | ||
| R2 Marg. | 0.180 | ||
| ICC | 0.9 | ||
| ELPD | -25.9 | ||
| ELPD s.e. | 6.1 | ||
| LOOIC | 51.7 | ||
| LOOIC s.e. | 12.2 | ||
| WAIC | 49.6 | ||
| RMSE | 0.21 | ||
10 Further analyses
In addition to comparing each of the symbiont types against the control group of no symbionts, it might be interesting to investigate whether there are any differences between the predation protection provided by crabs and shrimp, as well as whether having both crabs and shrimp symbionts is different to only a single symbiont type.
These contrasts can be explored via specific contrasts.
| SYMBIONT | Crab vs Shrimp | One vs Both | None vs Symbiont |
|---|---|---|---|
| none | 0 | 0 | 1 |
| crab | 1 | 1/2 | -1/3 |
| shrimp | -1 | 1/2 | -1/3 |
| both | 0 | -1 | -1/3 |
Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead. contrast odds.ratio lower.HPD upper.HPD
none / crabs 4.99 0.004819 47.6
none / shrimp 9.91 0.021648 88.6
none / both 27.17 0.118648 287.4
crabs / shrimp 2.02 0.000632 47.8
crabs / both 5.80 0.009004 183.7
shrimp / both 2.77 0.008164 45.6
Point estimate displayed: median
Results are back-transformed from the log odds ratio scale
HPD interval probability: 0.95 # via tidy_draws
mckeon_em <-
mckeon_brm4 |>
emmeans(~SYMBIONT, type = "link") |>
pairs() |>
tidy_draws() |>
mutate(across(starts_with("contrast"), exp)) |>
summarise_draws(median,
HDInterval::hdi,
Pl = ~ mean(.x < 1),
Pg = ~ mean(.x > 1),
ROPE = ~ mean(.x < 1.1 & .x > 0.9)
)Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.# A tibble: 6 × 7
variable median lower upper Pl Pg ROPE
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 contrast none - crabs 4.99 0.00482 47.6 0.185 0.815 0.0233
2 contrast none - shrimp 9.91 0.0216 88.6 0.0767 0.923 0.0171
3 contrast none - both 27.2 0.119 287. 0.0142 0.986 0.00333
4 contrast crabs - shrimp 2.02 0.000632 47.8 0.340 0.660 0.0429
5 contrast crabs - both 5.80 0.00900 184. 0.136 0.864 0.0288
6 contrast shrimp - both 2.77 0.00816 45.6 0.222 0.778 0.0408 # or via gather_emmeans_draws()
mckeon_em <-
mckeon_brm4 |>
emmeans(~SYMBIONT, type = "link") |>
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.# A tibble: 6 × 9
contrast y ymin ymax .width .point .interval Pl Pg
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <dbl> <dbl>
1 crabs - both 5.80 0.00900 183. 0.95 median hdci 0.136 0.864
2 crabs - shrimp 2.02 0.000632 47.4 0.95 median hdci 0.340 0.660
3 none - both 27.2 0.119 287. 0.95 median hdci 0.0142 0.986
4 none - crabs 4.99 0.00482 47.6 0.95 median hdci 0.185 0.815
5 none - shrimp 9.91 0.0216 88.5 0.95 median hdci 0.0767 0.923
6 shrimp - both 2.77 0.00816 45.6 0.95 median hdci 0.222 0.778## On a probability scale
mckeon_em <-
mckeon_brm4 |>
emmeans(~SYMBIONT, type = "link") |>
regrid() |>
pairs() |>
tidy_draws() |>
summarise_draws(median,
HDInterval::hdi,
Pl = ~ mean(.x < 0),
Pg = ~ mean(.x > 0)
)Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.# A tibble: 6 × 6
variable median lower upper Pl Pg
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 contrast none - crabs 0.139 -0.186 0.643 0.185 0.815
2 contrast none - shrimp 0.265 -0.0636 0.769 0.0767 0.923
3 contrast none - both 0.508 0.0567 0.898 0.0142 0.986
4 contrast crabs - shrimp 0.0958 -0.409 0.652 0.340 0.660
5 contrast crabs - both 0.309 -0.237 0.836 0.136 0.864
6 contrast shrimp - both 0.190 -0.289 0.762 0.222 0.778# or via gather_emmeans_draws()
## mckeon_em <- mckeon_brm4 |>
## emmeans(~SYMBIONT, type='link') |>
## pairs() |>
## gather_emmeans_draws() |>
## mutate(Eff=exp(.value),
## PEff=100*(Eff-1))#,
## # Prob = plogis(.value))
## mckeon_em |> head()
## mckeon_em |>
## group_by(contrast) |>
## dplyr::select(contrast, Eff) |>
## median_hdi()
## mckeon_em |>
## group_by(contrast) |>
## summarize(Prob=sum(Eff>1)/n())
## On a probability scale
mckeon_em <- mckeon_brm4 |>
emmeans(~SYMBIONT, type = "link") |>
regrid() |>
pairs() |>
gather_emmeans_draws() |>
mutate(Eff = .value) # ,Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.# A tibble: 6 × 6
# Groups: contrast [1]
contrast .chain .iteration .draw .value Eff
<chr> <int> <int> <int> <dbl> <dbl>
1 none - crabs NA NA 1 -0.0281 -0.0281
2 none - crabs NA NA 2 0.597 0.597
3 none - crabs NA NA 3 0.0459 0.0459
4 none - crabs NA NA 4 0.111 0.111
5 none - crabs NA NA 5 -0.0693 -0.0693
6 none - crabs NA NA 6 0.387 0.387 # A tibble: 6 × 7
contrast Eff .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 crabs - both 0.309 -0.215 0.859 0.95 median hdi
2 crabs - shrimp 0.0958 -0.409 0.652 0.95 median hdi
3 none - both 0.508 0.0492 0.891 0.95 median hdi
4 none - crabs 0.139 -0.142 0.689 0.95 median hdi
5 none - shrimp 0.265 -0.0636 0.769 0.95 median hdi
6 shrimp - both 0.190 -0.320 0.732 0.95 median hdi Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.# A tibble: 4 × 7
SYMBIONT P .lower .upper .width .point .interval
<fct> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 none 0.944 0.788 0.999 0.95 median hdci
2 crabs 0.794 0.275 1.000 0.95 median hdci
3 shrimp 0.663 0.189 0.999 0.95 median hdci
4 both 0.421 0.0241 0.881 0.95 median hdci cmat <- cbind(
"Crab vs shrimp" = c(0, 1, -1, 0),
"Both vs One" = c(0, -1 / 2, -1 / 2, 1),
"Any vs None" = c(-1, 1 / 3, 1 / 3, 1 / 3)
)
mckeon_em <- mckeon_brm4 |>
emmeans(~SYMBIONT, type = "response") |>
contrast(method = list(cmat))Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.mckeon_em <-
mckeon_brm4 |>
emmeans(~SYMBIONT, type = "link") |>
contrast(method = list(cmat)) |>
tidy_draws() |>
mutate(across(starts_with("contrast"), exp)) |>
summarise_draws(median,
HDInterval::hdi,
Pl = ~ mean(.x < 1),
Pg = ~ mean(.x > 1)
)Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.# A tibble: 3 × 6
variable median lower upper Pl Pg
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 contrast Crab vs shrimp 2.02 0.000632 47.8 0.340 0.660
2 contrast Both vs One 0.237 0.000305 2.15 0.861 0.139
3 contrast Any vs None 0.0942 0.00194 0.753 0.969 0.0308## or via gather_emmeans_draws
mckeon_em <-
mckeon_brm4 |>
emmeans(~SYMBIONT, type = "link") |>
contrast(method = list(cmat)) |>
gather_emmeans_draws() |>
mutate(Fit = exp(.value)) |>
summarise(median_hdci(Fit),
Pl = mean(Fit < 0),
Pg = mean(Fit > 0)
)Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead.# A tibble: 3 × 9
contrast y ymin ymax .width .point .interval Pl Pg
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <dbl> <dbl>
1 Any vs None 0.0942 0.00194 0.752 0.95 median hdci 0 1
2 Both vs One 0.237 0.000305 2.15 0.95 median hdci 0 1
3 Crab vs shrimp 2.02 0.000632 47.4 0.95 median hdci 0 1Warning: Method 'parse_bf' is deprecated. Please use 'brmsterms' instead. SYMBIONT response lower.HPD upper.HPD
none 0.9443812 0.7882106 0.9989135
crabs 0.7944272 0.2751096 0.9997865
shrimp 0.6632351 0.1895295 0.9995813
both 0.4212663 0.0091733 0.8548816
Point estimate displayed: median
Results are back-transformed from the logit scale
HPD interval probability: 0.95 ggplot(newdata, aes(y = response, x = SYMBIONT)) +
geom_pointrange(aes(ymin = lower.HPD, ymax = upper.HPD))
Estimate Est.Error Q2.5 Q97.5
R2 0.1763405 0.09054808 0.01394755 0.3286097 y ymin ymax .width .point .interval
1 0.1798171 0.008074149 0.321401 0.95 median hdci y ymin ymax .width .point .interval
1 0.4706625 0.1406366 0.734526 0.95 median hdci## for random intercept/slope model
mckeon_brm4 |>
bayes_R2(re.form = ~ (SYMBIONT | BLOCK), summary = FALSE) |>
median_hdci() y ymin ymax .width .point .interval
1 0.7080777 0.5762833 0.8229876 0.95 median hdci11 References
Mckeon, Seabird, Adrian Stier, Shelby Mcilroy, and Benjamin Bolker. 2012. “Multiple Defender Effects: Synergistic Coral Defense by Mutualist Crustaceans.” Oecologia 169 (February): 1095–1103. https://doi.org/10.1007/s00442-012-2275-2.