Estimating Distributional Models with
brms
Paul Bürkner
2024-03-19
Introduction
This vignette provides an introduction on how to fit distributional
regression models with brms. We use the term
distributional model to refer to a model, in which we can
specify predictor terms for all parameters of the assumed response
distribution. In the vast majority of regression model implementations,
only the location parameter (usually the mean) of the response
distribution depends on the predictors and corresponding regression
parameters. Other parameters (e.g., scale or shape parameters) are
estimated as auxiliary parameters assuming them to be constant across
observations. This assumption is so common that most researchers
applying regression models are often (in my experience) not aware of the
possibility of relaxing it. This is understandable insofar as relaxing
this assumption drastically increase model complexity and thus makes
models hard to fit. Fortunately, brms uses
Stan on the backend, which is an incredibly flexible
and powerful tool for estimating Bayesian models so that model
complexity is much less of an issue.
Suppose we have a normally distributed response variable. Then, in
basic linear regression, we specify a predictor term \(\eta_{\mu}\) for the mean parameter \(\mu\) of the normal distribution. The
second parameter of the normal distribution – the residual standard
deviation \(\sigma\) – is assumed to be
constant across observations. We estimate \(\sigma\) but do not try to predict
it. In a distributional model, however, we do exactly this by specifying
a predictor term \(\eta_{\sigma}\) for
\(\sigma\) in addition to the predictor
term \(\eta_{\mu}\). Ignoring
group-level effects for the moment, the linear predictor of a parameter
\(\theta\) for observation \(n\) has the form
\[\eta_{\theta n} = \sum_{i =
1}^{K_{\theta}} b_{\theta i} x_{\theta i n}\] where \(x_{\theta i n}\) denotes the value of the
\(i\)th predictor of parameter \(\theta\) for observation \(n\) and \(b_{\theta i}\) is the \(i\)th regression coefficient of parameter
\(\theta\). A distributional normal
model with response variable \(y\) can
then be written as
\[y_n \sim \mathcal{N}\left(\eta_{\mu n},
\, \exp(\eta_{\sigma n}) \right)\] We used the exponential
function around \(\eta_{\sigma}\) to
reflect that \(\sigma\) constitutes a
standard deviation and thus only takes on positive values, while a
linear predictor can be any real number.
A simple distributional model
Unequal variance models are possibly the most simple, but
nevertheless very important application of distributional models.
Suppose we have two groups of patients: One group receives a treatment
(e.g., an antidepressive drug) and another group receives placebo. Since
the treatment may not work equally well for all patients, the symptom
variance of the treatment group may be larger than the symptom variance
of the placebo group after some weeks of treatment. For simplicity,
assume that we only investigate the post-treatment values.
group <- rep(c("treat", "placebo"), each = 30)
symptom_post <- c(rnorm(30, mean = 1, sd = 2), rnorm(30, mean = 0, sd = 1))
dat1 <- data.frame(group, symptom_post)
head(dat1)
group symptom_post
1 treat -0.4659072
2 treat -2.9203450
3 treat -0.9915646
4 treat 1.2631294
5 treat 1.3677749
6 treat 4.3684730
The following model estimates the effect of group on
both the mean and the residual standard deviation of the normal response
distribution.
fit1 <- brm(bf(symptom_post ~ group, sigma ~ group),
data = dat1, family = gaussian())
Useful summary statistics and plots can be obtained via
summary(fit1)
plot(fit1, N = 2, ask = FALSE)
plot(conditional_effects(fit1), points = TRUE)
The population-level effect sigma_grouptreat, which is
the contrast of the two residual standard deviations on the log-scale,
reveals that the variances of both groups are indeed different. This
impression is confirmed when looking at the
conditional_effects of group. Going one step
further, we can compute the residual standard deviations on the original
scale using the hypothesis method.
hyp <- c("exp(sigma_Intercept) = 0",
"exp(sigma_Intercept + sigma_grouptreat) = 0")
hypothesis(fit1, hyp)
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (exp(sigma_Interc... = 0 1.42 0.19 1.11 1.87 NA NA *
2 (exp(sigma_Interc... = 0 1.99 0.27 1.56 2.61 NA NA *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
We may also directly compare them and plot the posterior distribution
of their difference.
hyp <- "exp(sigma_Intercept + sigma_grouptreat) > exp(sigma_Intercept)"
(hyp <- hypothesis(fit1, hyp))
Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (exp(sigma_Interc... > 0 0.57 0.33 0.04 1.14 26.21 0.96 *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
Indeed, the residual standard deviation of the treatment group seems
to larger than that of the placebo group. Moreover the magnitude of this
difference is pretty similar to what we expected due to the values we
put into the data simulations.
Zero-Inflated Models
Another important application of the distributional regression
framework are so called zero-inflated models. These models are helpful
whenever there are more zeros in the response variable than one would
naturally expect. For example, if one seeks to predict the number of
cigarettes people smoke per day and also includes non-smokers, there
will be a huge amount of zeros which, when not modeled appropriately,
can seriously distort parameter estimates. Here, we consider an example
dealing with the number of fish caught by various groups of people. On
the UCLA website (), the data are described as follows: “The state
wildlife biologists want to model how many fish are being caught by
fishermen at a state park. Visitors are asked how long they stayed, how
many people were in the group, were there children in the group and how
many fish were caught. Some visitors do not fish, but there is no data
on whether a person fished or not. Some visitors who did fish did not
catch any fish so there are excess zeros in the data because of the
people that did not fish.”
zinb <- read.csv("https://paul-buerkner.github.io/data/fish.csv")
head(zinb)
nofish livebait camper persons child xb zg count
1 1 0 0 1 0 -0.8963146 3.0504048 0
2 0 1 1 1 0 -0.5583450 1.7461489 0
3 0 1 0 1 0 -0.4017310 0.2799389 0
4 0 1 1 2 1 -0.9562981 -0.6015257 0
5 0 1 0 1 0 0.4368910 0.5277091 1
6 0 1 1 4 2 1.3944855 -0.7075348 0
As predictors we choose the number of people per group, the number of
children, as well as whether the group consists of campers. Many groups
may not even try catching any fish at all (thus leading to many zero
responses) and so we fit a zero-inflated Poisson model to the data. For
now, we assume a constant zero-inflation probability across
observations.
fit_zinb1 <- brm(count ~ persons + child + camper,
data = zinb, family = zero_inflated_poisson())
Again, we summarize the results using the usual methods.
Family: zero_inflated_poisson
Links: mu = log; zi = identity
Formula: count ~ persons + child + camper
Data: zinb (Number of observations: 250)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -1.01 0.17 -1.36 -0.68 1.00 2757 2751
persons 0.87 0.04 0.79 0.97 1.00 2921 2550
child -1.37 0.09 -1.56 -1.19 1.00 2746 2597
camper 0.80 0.09 0.62 0.98 1.00 3230 2471
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
zi 0.41 0.05 0.31 0.49 1.00 3023 2672
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).
plot(conditional_effects(fit_zinb1), ask = FALSE)
According to the parameter estimates, larger groups catch more fish,
campers catch more fish than non-campers, and groups with more children
catch less fish. The zero-inflation probability zi is
pretty large with a mean of 41%. Please note that the probability of
catching no fish is actually higher than 41%, but parts of this
probability are already modeled by the Poisson distribution itself
(hence the name zero-inflation). If you want to treat all zeros
as originating from a separate process, you can use hurdle models
instead (not shown here).
Now, we try to additionally predict the zero-inflation probability by
the number of children. The underlying reasoning is that we expect
groups with more children to not even try catching fish. Most children
are just terribly bad at waiting for hours until something happens. From
a purely statistical perspective, zero-inflated (and hurdle)
distributions are a mixture of two processes and predicting both parts
of the model is natural and often very reasonable to make full use of
the data.
fit_zinb2 <- brm(bf(count ~ persons + child + camper, zi ~ child),
data = zinb, family = zero_inflated_poisson())
Family: zero_inflated_poisson
Links: mu = log; zi = logit
Formula: count ~ persons + child + camper
zi ~ child
Data: zinb (Number of observations: 250)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -1.09 0.18 -1.45 -0.74 1.00 2825 2524
zi_Intercept -0.97 0.27 -1.53 -0.47 1.00 3480 2175
persons 0.89 0.05 0.81 0.99 1.00 2979 2559
child -1.18 0.09 -1.36 -1.00 1.00 3210 2770
camper 0.78 0.09 0.60 0.96 1.00 3351 2765
zi_child 1.22 0.27 0.69 1.80 1.00 3779 3126
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).
plot(conditional_effects(fit_zinb2), ask = FALSE)
To transform the linear predictor of zi into a
probability, brms applies the logit-link:
\[logit(zi) =
\log\left(\frac{zi}{1-zi}\right) = \eta_{zi}\]
The logit-link takes values within \([0,
1]\) and returns values on the real line. Thus, it allows the
transition between probabilities and linear predictors.
According to the model, trying to fish with children not only
decreases the overall number fish caught (as implied by the Poisson part
of the model) but also drastically increases your change of catching no
fish at all (as implied by the zero-inflation part) most likely because
groups with more children are not even trying.
Additive Distributional Models
In the examples so far, we did not have multilevel data and thus did
not fully use the capabilities of the distributional regression
framework of brms. In the example presented below, we
will not only show how to deal with multilevel data in distributional
models, but also how to incorporate smooth terms (i.e., splines) into
the model. In many applications, we have no or only a very vague idea
how the relationship between a predictor and the response looks like. A
very flexible approach to tackle this problems is to use splines and let
them figure out the form of the relationship. For illustration purposes,
we simulate some data with the mgcv package, which is
also used in brms to prepare smooth terms.
dat_smooth <- mgcv::gamSim(eg = 6, n = 200, scale = 2, verbose = FALSE)
Gu & Wahba 4 term additive model
y x0 x1 x2 x3 f
1 16.654758 0.8283741 0.76433527 0.29431102 0.4183742 16.482587
2 15.117487 0.1894057 0.59874698 0.07591961 0.2629901 12.419708
3 14.693045 0.8315630 0.83349578 0.80335710 0.1967228 16.346423
4 22.152958 0.3372437 0.13709704 0.40472855 0.2868083 19.188167
5 12.926583 0.4825321 0.32762949 0.34206137 0.2021112 13.128004
6 4.346698 0.1004876 0.02144073 0.75854505 0.1724274 9.563716
The data contains the predictors x0 to x3
as well as the grouping factor fac indicating the nested
structure of the data. We predict the response variable y
using smooth terms of x1 and x2 and a varying
intercept of fac. In addition, we assume the residual
standard deviation sigma to vary by a smoothing term of
x0 and a varying intercept of fac.
fit_smooth1 <- brm(
bf(y ~ s(x1) + s(x2) + (1|fac), sigma ~ s(x0) + (1|fac)),
data = dat_smooth, family = gaussian(),
chains = 2, control = list(adapt_delta = 0.95)
)
Family: gaussian
Links: mu = identity; sigma = log
Formula: y ~ s(x1) + s(x2) + (1 | fac)
sigma ~ s(x0) + (1 | fac)
Data: dat_smooth (Number of observations: 200)
Draws: 2 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 2000
Smoothing Spline Hyperparameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sds(sx1_1) 2.80 1.82 0.60 7.33 1.00 1179 1180
sds(sx2_1) 19.52 5.27 11.62 32.05 1.00 951 1116
sds(sigma_sx0_1) 1.21 1.27 0.02 4.84 1.00 525 787
Multilevel Hyperparameters:
~fac (Number of levels: 4)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 4.95 2.10 2.34 10.85 1.00 995 1227
sd(sigma_Intercept) 0.21 0.26 0.01 0.90 1.00 417 593
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 15.04 2.29 10.32 19.63 1.00 757 915
sigma_Intercept 0.69 0.17 0.36 1.03 1.00 828 578
sx1_1 12.80 5.58 1.49 24.14 1.00 1265 1175
sx2_1 39.86 13.93 12.20 67.25 1.00 1297 1466
sigma_sx0_1 -0.50 2.71 -6.51 5.52 1.00 1108 694
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).
plot(conditional_effects(fit_smooth1), points = TRUE, ask = FALSE)
This model is likely an overkill for the data at hand, but nicely
demonstrates the ease with which one can specify complex models with
brms and to fit them using Stan on the
backend.