# NOT RUN {
# Poisson regression for the number of seizures in epileptic patients
# using normal priors for population-level effects
# and half-cauchy priors for standard deviations of group-level effects
prior1 <- prior(normal(0,10), class = b) +
prior(cauchy(0,2), class = sd)
fit1 <- brm(count ~ zAge + zBase * Trt + (1|patient),
data = epilepsy, family = poisson(), prior = prior1)
# generate a summary of the results
summary(fit1)
# plot the MCMC chains as well as the posterior distributions
plot(fit1, ask = FALSE)
# predict responses based on the fitted model
head(predict(fit1))
# plot conditional effects for each predictor
plot(conditional_effects(fit1), ask = FALSE)
# investigate model fit
loo(fit1)
pp_check(fit1)
# Ordinal regression modeling patient's rating of inhaler instructions
# category specific effects are estimated for variable 'treat'
fit2 <- brm(rating ~ period + carry + cs(treat),
data = inhaler, family = sratio("logit"),
prior = set_prior("normal(0,5)"), chains = 2)
summary(fit2)
plot(fit2, ask = FALSE)
WAIC(fit2)
# Survival regression modeling the time between the first
# and second recurrence of an infection in kidney patients.
fit3 <- brm(time | cens(censored) ~ age * sex + disease + (1|patient),
data = kidney, family = lognormal())
summary(fit3)
plot(fit3, ask = FALSE)
plot(conditional_effects(fit3), ask = FALSE)
# Probit regression using the binomial family
ntrials <- sample(1:10, 100, TRUE)
success <- rbinom(100, size = ntrials, prob = 0.4)
x <- rnorm(100)
data4 <- data.frame(ntrials, success, x)
fit4 <- brm(success | trials(ntrials) ~ x, data = data4,
family = binomial("probit"))
summary(fit4)
# Non-linear Gaussian model
fit5 <- brm(
bf(cum ~ ult * (1 - exp(-(dev/theta)^omega)),
ult ~ 1 + (1|AY), omega ~ 1, theta ~ 1,
nl = TRUE),
data = loss, family = gaussian(),
prior = c(
prior(normal(5000, 1000), nlpar = "ult"),
prior(normal(1, 2), nlpar = "omega"),
prior(normal(45, 10), nlpar = "theta")
),
control = list(adapt_delta = 0.9)
)
summary(fit5)
conditional_effects(fit5)
# Normal model with heterogeneous variances
data_het <- data.frame(
y = c(rnorm(50), rnorm(50, 1, 2)),
x = factor(rep(c("a", "b"), each = 50))
)
fit6 <- brm(bf(y ~ x, sigma ~ 0 + x), data = data_het)
summary(fit6)
plot(fit6)
conditional_effects(fit6)
# extract estimated residual SDs of both groups
sigmas <- exp(as.data.frame(fit6, variable = "^b_sigma_", regex = TRUE))
ggplot(stack(sigmas), aes(values)) +
geom_density(aes(fill = ind))
# Quantile regression predicting the 25%-quantile
fit7 <- brm(bf(y ~ x, quantile = 0.25), data = data_het,
family = asym_laplace())
summary(fit7)
conditional_effects(fit7)
# use the future package for more flexible parallelization
library(future)
plan(multiprocess)
fit7 <- update(fit7, future = TRUE)
# fit a model manually via rstan
scode <- make_stancode(count ~ Trt, data = epilepsy)
sdata <- make_standata(count ~ Trt, data = epilepsy)
stanfit <- rstan::stan(model_code = scode, data = sdata)
# feed the Stan model back into brms
fit8 <- brm(count ~ Trt, data = epilepsy, empty = TRUE)
fit8$fit <- stanfit
fit8 <- rename_pars(fit8)
summary(fit8)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab