library(tidyverse)
library(brms)
library(posterior)
library(bayesplot)
library(ggdist)
library(patchwork)
color_scheme_set(scheme = "red")
theme_set(theme_ggdist())
prior_pred_lines <- function(post) {
map(.x = seq_len(nrow(post)),
.f = function(ii, post){
p_Intercept <- post$b_Intercept[ii]
p_shape <- post$shape[ii]
tibble(ii = ii,
x = seq(0, 30, length.out = 100),
y = dgamma(x, shape = exp(p_Intercept), rate = 1 / p_shape))
}, post = post) |>
list_rbind()
}Gamma
Resources
- https://discourse.mc-stan.org/t/output-for-family-in-brms/9379
- https://discourse.mc-stan.org/t/decoding-smooth-fits-in-brms/34473
- https://stats.stackexchange.com/questions/96972/how-to-interpret-parameters-in-glm-with-family-gamma
Gamma distribution
Shape (\(\theta\)) and Scale (\(k\))
\[PDF(x)=\frac{1}{\Gamma(k)~\theta^k} x^{k - 1}~e^{-x/\theta}\]
Shape (\(\alpha = \theta\)) and Rate (\(\beta\))
\[PDF(x) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha - 1}~e^{-\beta x}\]
Inverse scale parameter = rate
\[\beta = \frac{1}{\theta} = \frac{1}{\alpha}\]
\[\bar{x} = \mu = \frac{\alpha}{\beta}\]
Data
set.seed(3475934)
x_mean <- 4
alpha <- 12
(beta <- alpha / x_mean)[1] 3n <- 1e4
GD <- tibble(x = rgamma(n, shape = alpha, rate = beta))
alpha / beta[1] 4mean(GD$x)[1] 4.010262\[Var(x) = \frac{\alpha}{\beta^2}\]
alpha / (beta^2)[1] 1.333333var(GD$x)[1] 1.327161\[Var(x) = \phi~\mu^2\]
where \(\phi\) is a dispersion parameter
(alpha / (beta^2)) / ((alpha / beta)^2)[1] 0.08333333GLM
fm_glm <- glm(x ~ 1, data = GD, family = Gamma(link = log))
summary(fm_glm)
Call:
glm(formula = x ~ 1, family = Gamma(link = log), data = GD)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.388857 0.002873 483.5 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.08252358)
Null deviance: 842.49 on 9999 degrees of freedom
Residual deviance: 842.49 on 9999 degrees of freedom
AIC: 30718
Number of Fisher Scoring iterations: 4Mean
\[\mu = \exp(\alpha + \beta~x)\]
(fm_beta <- coef(fm_glm))(Intercept)
1.388857 exp(fm_beta)(Intercept)
4.010262 alpha / beta[1] 4# alpha
1 / MASS::gamma.dispersion(fm_glm)[1] 12.03391Priors
get_prior(x ~ 1,
family = Gamma(link = log),
data = GD) prior class coef group resp dpar nlpar lb ub source
student_t(3, 1.4, 2.5) Intercept default
gamma(0.01, 0.01) shape 0 defaultPrior prediction - default
priors <- c(prior(student_t(3, 1.4, 2.5), class = "Intercept"),
prior(gamma(0.01, 0.01), class = "shape"))
PP <- brm(x ~ 1,
family = Gamma(link = log),
prior = priors,
data = GD,
sample_prior = "only",
backend = "cmdstanr",
refresh = 0)Running MCMC with 4 sequential chains...
Chain 1 finished in 0.0 seconds.
Chain 2 finished in 0.0 seconds.
Chain 3 finished in 0.0 seconds.
Chain 4 finished in 0.1 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.0 seconds.
Total execution time: 0.8 seconds.(post <- as_draws_df(PP) |>
slice_sample(n = 200))# A draws_df: 180 iterations, 4 chains, and 5 variables
b_Intercept shape Intercept lprior lp__
1 1.87 6.1e+00 1.87 -8.4 -6.6
2 4.32 2.0e-207 4.32 463.9 -12.1
3 0.81 1.8e-72 0.81 156.9 -8.3
4 5.61 2.9e-08 5.61 9.3 -8.1
5 18.50 2.1e-21 18.50 34.9 -12.7
6 2.33 6.2e-251 2.33 563.7 -12.4
7 5.45 9.1e-105 5.45 229.4 -10.2
8 7.34 1.8e-44 7.34 91.1 -9.7
9 -2.43 1.9e-19 -2.43 35.0 -8.1
10 6.93 2.0e-16 6.93 27.3 -8.9
# ... with 190 more draws
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}ggplot() +
geom_line(data = prior_pred_lines(post), aes(x, y, group = ii),
color = "firebrick4", alpha = 0.5) +
geom_dotsinterval(data = GD, aes(x))
Prior prediction - improved
(pr_mean <- log(mean(GD$x))) |> round(2)[1] 1.39priors <- c(prior(student_t(3, 1.4, 2.5), class = "Intercept"),
prior(gamma(1.4, 1/1.4), class = "shape"))
PP <- brm(x ~ 1,
family = Gamma(link = log),
prior = priors,
data = GD,
backend = "cmdstanr",
sample_prior = "only",
refresh = 0)Start samplingRunning MCMC with 4 sequential chains...
Chain 1 finished in 0.0 seconds.
Chain 2 finished in 0.0 seconds.
Chain 3 finished in 0.0 seconds.
Chain 4 finished in 0.0 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.0 seconds.
Total execution time: 0.6 seconds.(post <- as_draws_df(PP) |>
slice_sample(n = 200))# A draws_df: 188 iterations, 4 chains, and 5 variables
b_Intercept shape Intercept lprior lp__
1 17.42 1.71 17.42 -8.6 -8.1
2 2.63 1.17 2.63 -3.2 -3.0
3 0.50 6.27 0.50 -6.1 -4.3
4 6.25 0.71 6.25 -4.5 -4.9
5 -5.96 1.20 -5.96 -5.8 -5.6
6 0.44 4.31 0.44 -4.9 -3.4
7 3.91 0.94 3.91 -3.5 -3.6
8 10.46 3.19 10.46 -7.4 -6.3
9 2.83 1.23 2.83 -3.3 -3.1
10 -2.25 4.07 -2.25 -5.7 -4.3
# ... with 190 more draws
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}ggplot() +
geom_line(data = prior_pred_lines(post), aes(x, y, group = ii),
color = "firebrick4", alpha = 0.5) +
geom_dotsinterval(data = GD, aes(x))
brms model
fm <- brm(x ~ 1,
family = Gamma(link = log),
data = GD,
prior = priors,
backend = "cmdstanr",
cores = 4,
init = list(list(Intercept = 1, shape = 2),
list(Intercept = 1, shape = 2),
list(Intercept = 1, shape = 2),
list(Intercept = 1, shape = 2)),
control = list(step_size = 0.01),
refresh = 0)Start samplingRunning MCMC with 4 parallel chains...
Chain 3 finished in 6.7 seconds.
Chain 4 finished in 6.8 seconds.
Chain 1 finished in 7.0 seconds.
Chain 2 finished in 7.1 seconds.
All 4 chains finished successfully.
Mean chain execution time: 6.9 seconds.
Total execution time: 7.2 seconds.prior_summary(fm) prior class coef group resp dpar nlpar lb ub source
student_t(3, 1.4, 2.5) Intercept user
gamma(1.4, 1/1.4) shape 0 usermcmc_combo(fm, pars = c("b_Intercept", "shape"),
combo = c("dens_overlay", "rank_overlay"))
pp_check(fm, ndraws = 100)
ps <- posterior_summary(fm)(mean_est <- exp(ps[1, 1]))[1] 4.01048mean(GD$x)[1] 4.010262