# NOT RUN {
## demonstrate how to fit a beta-binomial model
## generate some fake data
phi <- 0.7
n <- 300
z <- rnorm(n, sd = 0.2)
ntrials <- sample(1:10, n, replace = TRUE)
eta <- 1 + z
mu <- exp(eta) / (1 + exp(eta))
a <- mu * phi
b <- (1 - mu) * phi
p <- rbeta(n, a, b)
y <- rbinom(n, ntrials, p)
dat <- data.frame(y, z, ntrials)
# define a custom family
beta_binomial2 <- custom_family(
"beta_binomial2", dpars = c("mu", "phi"),
links = c("logit", "log"), lb = c(NA, 0),
type = "int", vars = "trials[n]"
)
# define the corresponding Stan density function
stan_funs <- "
real beta_binomial2_lpmf(int y, real mu, real phi, int N) {
return beta_binomial_lpmf(y | N, mu * phi, (1 - mu) * phi);
}
"
# fit the model
fit <- brm(y | trials(ntrials) ~ z, data = dat,
family = beta_binomial2, stan_funs = stan_funs)
summary(fit)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab