library(tidyverse)
library(brms)
library(tidybayes)
library(marginaleffects)
# Example data with proportions
example_data <- tribble(
~answer, ~n, ~total,
"Very familiar", 88, 205,
"Somewhat familiar", 80, 249,
"Not familiar at all", 19, 85
)
example_data %>%
mutate(prop = n / total)
#> # A tibble: 3 × 4
#> answer n total prop
#> <chr> <dbl> <dbl> <dbl>
#> 1 Very familiar 88 205 0.429
#> 2 Somewhat familiar 80 249 0.321
#> 3 Not familiar at all 19 85 0.224
# Model with identity-linked binomial family
mod <- brm(
bf(n | trials(total) ~ 0 + answer),
data = example_data,
family = binomial(link = "identity"),
prior = c(prior(beta(2, 2), class = b, ub = 1, lb = 0)),
backend = "cmdstanr"
)
# Same proportions!
mod
#> Family: binomial
#> Links: mu = identity
#> Formula: n | trials(total) ~ 0 + answer
#> Data: example_data (Number of observations: 3)
#> Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#> total post-warmup draws = 4000
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> answerNotfamiliaratall 0.24 0.04 0.16 0.33 1.00 4377
#> answerSomewhatfamiliar 0.32 0.03 0.27 0.38 1.00 4624
#> answerVeryfamiliar 0.43 0.03 0.37 0.50 1.00 4704
#>
#> Draws were sampled using sample(hmc). 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).
# Pairwise differences!
mod %>%
gather_draws(`^b_answer.*`, regex = TRUE) %>%
ungroup() %>%
compare_levels(.value, by = .variable) %>%
ggplot(aes(x = .value, fill = .variable)) +
stat_halfeye() +
guides(fill = "none") +
facet_wrap(vars(.variable), ncol = 1)
# Pairwise differences with marginaleffects
mod %>%
avg_comparisons(variables = list(answer = "pairwise"),
type = "link") %>%
posterior_draws() %>%
ggplot(aes(x = draw, fill = contrast)) +
stat_halfeye() +
guides(fill = "none") +
facet_wrap(vars(contrast), ncol = 1)
Created on 2023-05-09 with reprex v2.0.2