# source("renv/activate.R")
# Sys.setenv(RENV_DOWNLOAD_METHOD = "curl")
# options(renv.install.staged = FALSE)
# renv::settings$use.cache(TRUE)
# Sys.setenv(RENV_INSTALL_PAK_ENABLED = "TRUE")Practical Bayesian Modeling with Stan and brms
1 Docker (optional)
Run this on your terminal to start an RStudio Server with all the necessary packages installed.
# run in afec-bayes-2025 directory
docker run -d -p 8787:8787 \
-e PASSWORD=123456 \
-v $(pwd):/home/rstudio/workspace \
mattocci/afec-bayes:prodAccess 127.0.0.1:8787 and log in with the username rstudio and the password 123456 (or your own password).
Note: Because you’re running this container on your own machine and mapping only to 127.0.0.1:8787, the RStudio Server isn’t accessible from outside your computer. So, using a simple password like 123456 is fine and safe for local development.
2 renv settings (Optional)
Comment out the following lines in .Rprofile if you are not using renv.
3 Setup
library(brms)
library(cmdstanr)
library(tidyverse)
library(bayesplot)
library(patchwork)
cmdstanr::cmdstan_path()[1] "/home/rstudio/.cmdstan/cmdstan-2.37.0"my_theme <- theme_bw() + theme(
axis.text = element_text(size = 14),
strip.text = element_text(size = 14)
)
theme_set(my_theme)4 Normal model
Leaf mass per area (LMA) is an important trait that reflects plant strategies. LMA for interspecific data often shows a log-normal distribution.
\[ \text{log}(y_i) \sim N(\mu, \sigma) \]
Q: What are the mean (\(\mu\)) and standard deviation (\(\sigma\)) of log LMA?
4.1 Dummy data
set.seed(123)
n <- 100
mu <- 4.6
sigma <- 0.7
y <- rnorm(n, mean = mu, sd = sigma)
hist(y)
cmdstan uses list format for data input.
normal_list <- list(
y = y,
N = length(y)
)brms uses dataframe (tibble) format for data input like lme4::lmer and stats::glm.
normal_df <- tibble(
y = y,
N = length(y)
)4.2 CmdStan
We need to write and save Stan programs (.stan files).
stan/normal.stan
// observed data
data {
int<lower=0> N;
vector[N] y;
}
// parameters to be estimated
parameters {
real mu;
real<lower=0> sigma;
}
model {
// likelihood
y ~ normal(mu, sigma);
// priors
mu ~ normal(0, 10);
sigma ~ cauchy(0, 2.5);
} stan/normal_vague.stan
data {
int<lower=0> N;
vector[N] y;
}
parameters {
real mu;
real<lower=0> sigma;
}
model {
y ~ normal(mu, sigma);
mu ~ normal(0, 1e+4); // Vague prior for mu
sigma ~ normal(0, 1e+4); // Vague prior for sigma
} These commands compile the models, translating Stan code to C++ and then into machine-executable files. This takes a while (about 20 secs to 1 min) the first time only.
normal_mod <- cmdstan_model("stan/normal.stan")
normal_vague_mod <- cmdstan_model("stan/normal_vague.stan")# Persist CSV outputs to avoid /tmp cleanup during render
normal_out_dir <- file.path("assets", "cmdstan", "normal")
if (!dir.exists(normal_out_dir)) dir.create(normal_out_dir, recursive = TRUE)
normal_fit <- normal_mod$sample(
data = normal_list,
seed = 123,
chains = 4,
parallel_chains = 4,
iter_warmup = 1000, # number of warmup iterations
iter_sampling = 1000, # number of sampling iterations
refresh = 500, # print update every 500 iters
output_dir = normal_out_dir
)Running MCMC with 4 parallel chains...
Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 1 Iteration: 2000 / 2000 [100%] (Sampling) Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp10tuT6/model-23c7595ff71.stan', line 15, column 2 to column 24)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
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.3 seconds.We don’t have to recompile the model when we want to change the number of iterations, chains, or data.
normal_vague_fit <- normal_vague_mod$sample(
data = normal_list,
seed = 123,
chains = 4,
parallel_chains = 4,
refresh = 500 # print update every 500 iters
)4.2.1 Summary
We use this output when writing manuscripts and for model diagnostics.
normal_fit$summary()# A tibble: 3 × 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ -6.35 -6.04 0.997 0.731 -8.36 -5.38 1.00 1801. 2192.
2 mu 4.66 4.66 0.0657 0.0661 4.56 4.77 1.00 3300. 2446.
3 sigma 0.648 0.646 0.0466 0.0449 0.577 0.730 1.00 3181. 2516.lp__: log posterior.- \(\text{log}\; p(\theta \mid \text{data}) \propto \text{log}\; p(\text{data} \mid \theta) + \text{log}\; p(\theta)\)
rhat(Gelman-Rubin statistic): should be close to 1.- Measures the convergence of the chains.
rhat> 1.1: Definitely bad.- 1.01 <
rhat< 1.1: Suspicious. rhat<= 1.01: Good.
ess_bulkandess_tail: Effective sample size for bulk and tail of the posterior distribution.ess_bulk: sampling efficiency for the bulk posteriors (e.g., mean, SD).ess_tail: sampling efficiency for the tails (e.g., quantiles like 2.5%, 97.5%).ess_bulk,ess_tail> 400: Good.
4.2.2 Draws (posterior samples)
For each parameter, we have 1000 iterations \(\times\) 4 chains = 4000 posteriors. We use this for visualization and diagnostics.
normal_fit$draws()# A draws_array: 1000 iterations, 4 chains, and 3 variables
, , variable = lp__
chain
iteration 1 2 3 4
1 -7.9 -5.5 -6.7 -6.2
2 -7.6 -5.4 -8.9 -6.8
3 -6.7 -5.3 -6.4 -8.3
4 -6.1 -5.6 -5.9 -6.7
5 -7.8 -5.6 -6.3 -7.3
, , variable = mu
chain
iteration 1 2 3 4
1 4.6 4.6 4.6 4.7
2 4.6 4.7 4.8 4.7
3 4.6 4.7 4.7 4.8
4 4.6 4.7 4.6 4.7
5 4.8 4.7 4.7 4.7
, , variable = sigma
chain
iteration 1 2 3 4
1 0.75 0.63 0.72 0.65
2 0.57 0.65 0.59 0.57
3 0.69 0.64 0.61 0.57
4 0.64 0.62 0.66 0.58
5 0.65 0.62 0.70 0.58
# ... with 995 more iterationsTrace plots for diagnosing convergence and mixing of the chains.
normal_draws <- as_draws_df(normal_fit$draws())
color_scheme_set("viridis")
mcmc_trace(normal_draws, pars = c("mu", "sigma"))
Histograms of the posterior distributions of the parameters.
mcmc_hist(normal_draws, pars = c("mu", "sigma"))`stat_bin()` using `bins = 30`. Pick better value `binwidth`.
4.2.3 Diagnostic summary
normal_fit$diagnostic_summary()$num_divergent
[1] 0 0 0 0
$num_max_treedepth
[1] 0 0 0 0
$ebfmi
[1] 1.151656 1.101592 1.137121 1.100724num_divergent: indicates the number of iterations (sampling transitions) where the Hamiltonian trajectory failed (divergent transition).- Even 1-2 divergent transitions suggest that model may not be reliable, especially in hierarchical models.
- We need reparameterization or increase
adapt_deltato make the sampler take smaller steps.
num_max_treedepth: indicates the number of iterations where there are not enough leapfrog steps.- This can indicate that the model is complex or that the priors are too tight.
- We can increase
max_treedepthto allow more steps, but this can increase the computation time.
ebfmi: Energy-Bayesian Fraction of Missing Information- Measures whether the resampled momentum generates enough variation in energy to explore the posterior efficiently.
- Low
ebfmimeans the sampler may be stuck in tight regions and exploring poorly, even if Rhat and ESS look fine. - Guidelines:
ebfmi< 0.3: Bad- 0.3 <
ebfmi<= 1: Acceptable ebfmi>= 1: Good
4.2.4 Methods text example
“We estimated posterior distributions of all parameters using the Hamiltonian Monte Carlo (HMC) algorithm implemented in Stan (Carpenter et al., 2017) with weakly informative priors (Gelman et al., 2008). The HMC algorithm uses gradient information to propose new states in the Markov chain, leading to a more efficient exploration of the target distribution than traditional Markov chain Monte Carlo (MCMC) methods that rely on random proposals (Carpenter et al., 2017). This efficiency allows us to achieve convergence with fewer iterations than traditional MCMC methods. Four independent chains were run for 2000 iterations for each model with a warm-up of 1000 iterations. Convergence of the posterior distribution was assessed using the Gelman–Rubin statistic with a convergence threshold of 1.1 (Gelman et al., 2013), ensuring effective sample sizes greater than 400 (Vehtari et al., 2021), and by monitoring divergent transitions (Betancourt, 2016) for all parameters.”
4.3 brms
brms uses an lme4-like syntax.
normal_fit_brm0 <- brm(y ~ 1, family = gaussian(), data = normal_df)We can also specify priors for the parameters. The prior argument can be a bit tricky, since it is not always obvious which parameter names correspond to which parts of the model, especially when the model becomes more complex.
We can also control the number of warmup iterations, total iterations, and other sampling settings.
The cmdstanr backend is generally faster than the default rstan backend.
normal_fit_brm <- brm(
y ~ 1,
family = gaussian(),
data = normal_df,
prior = c(
prior(normal(0, 5), class = "Intercept"),
prior(cauchy(0, 2.5), class = "sigma")
),
seed = 123,
iter = 2000, # total iterations (warmup + post-warmup)
warmup = 1000,
chains = 4, cores = 4,
backend = "cmdstanr"
)In file included from stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/hypergeometric_pFq.hpp:11,
from stan/lib/stan_math/stan/math/prim/fun/hypergeometric_pFq.hpp:8,
from stan/lib/stan_math/stan/math/prim/fun/hypergeometric_2F1.hpp:19,
from stan/lib/stan_math/stan/math/prim/fun/grad_2F1.hpp:14,
from stan/lib/stan_math/stan/math/prim/fun.hpp:117,
from stan/lib/stan_math/stan/math/prim/constraint.hpp:4,
from stan/lib/stan_math/stan/math/rev/constraint.hpp:24,
from stan/lib/stan_math/stan/math/rev.hpp:11,
from stan/lib/stan_math/stan/math.hpp:19,
from stan/src/stan/model/model_header.hpp:6:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp: In instantiation of ‘std::pair<_ForwardIterator, _ForwardIterator> boost::math::detail::hypergeometric_pFq_checked_series_impl(const Seq&, const Seq&, const Real&, const Policy&, const Terminal&, long long int&) [with Seq = std::vector<double>; Real = long double; Policy = boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>; Terminal = iteration_terminator]’:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/hypergeometric_pFq.hpp:61:107: required from ‘typename boost::math::tools::promote_args<Real, typename Seq::value_type>::type boost::math::hypergeometric_pFq(const Seq&, const Seq&, const Real&, Real*, const Policy&) [with Seq = std::vector<double>; Real = double; Policy = policies::policy<policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy>; typename tools::promote_args<Real, typename Seq::value_type>::type = double; typename Seq::value_type = double]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/hypergeometric_pFq.hpp:72:35: required from ‘typename boost::math::tools::promote_args<Real, typename Seq::value_type>::type boost::math::hypergeometric_pFq(const Seq&, const Seq&, const Real&, Real*) [with Seq = std::vector<double>; Real = double; typename tools::promote_args<Real, typename Seq::value_type>::type = double; typename Seq::value_type = double]’
stan/lib/stan_math/stan/math/prim/fun/hypergeometric_pFq.hpp:53:41: required from ‘stan::return_type_t<T, L, U> stan::math::hypergeometric_pFq(const Ta&, const Tb&, const Tz&) [with Ta = Eigen::Matrix<double, -1, 1>; Tb = Eigen::Matrix<double, -1, 1>; Tz = double; stan::require_all_vector_st<std::is_arithmetic, Ta, Tb>* <anonymous> = 0; stan::require_arithmetic_t<Tz>* <anonymous> = 0; stan::return_type_t<T, L, U> = double]’
stan/lib/stan_math/stan/math/prim/fun/hypergeometric_3F2.hpp:130:28: required from ‘auto stan::math::hypergeometric_3F2(const Ta&, const Tb&, const Tz&) [with Ta = std::vector<double>; Tb = std::vector<double>; Tz = double; stan::require_all_vector_t<T_y, T_high>* <anonymous> = 0; stan::require_stan_scalar_t<T2>* <anonymous> = 0]’
stan/lib/stan_math/stan/math/prim/fun/hypergeometric_3F2.hpp:154:28: required from ‘auto stan::math::hypergeometric_3F2(const std::initializer_list<_Tp>&, const std::initializer_list<_Value>&, const Tz&) [with Ta = double; Tb = double; Tz = double; stan::require_all_stan_scalar_t<T, L, U>* <anonymous> = 0]’
stan/lib/stan_math/stan/math/rev/fun/inv_inc_beta.hpp:80:37: required from here
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp:122:28: note: parameter passing for argument of type ‘std::pair<long double, long double>’ when C++17 is enabled changed to match C++14 in GCC 10.1
122 | std::pair<Real, Real> hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, const Terminal& termination, long long& log_scale)
| ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~In file included from stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/factorials.hpp:17,
from stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/binomial.hpp:16,
from stan/lib/stan_math/stan/math/prim/fun/choose.hpp:7,
from stan/lib/stan_math/stan/math/prim/functor/hcubature.hpp:21,
from stan/lib/stan_math/stan/math/prim/functor.hpp:14,
from stan/lib/stan_math/stan/math/rev/fun/to_arena.hpp:7,
from stan/lib/stan_math/stan/math/rev/core/operator_division.hpp:15,
from stan/lib/stan_math/stan/math/rev/core/operator_divide_equal.hpp:5,
from stan/lib/stan_math/stan/math/rev/core.hpp:29,
from stan/src/stan/model/model_base.hpp:8,
from stan/src/stan/model/model_header.hpp:4:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/gamma.hpp: In instantiation of ‘boost::math::detail::upper_incomplete_gamma_fract<T>::result_type boost::math::detail::upper_incomplete_gamma_fract<T>::operator()() [with T = double; result_type = std::pair<double, double>]’:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/tools/fraction.hpp:259:20: required from ‘typename boost::math::tools::detail::fraction_traits<Gen>::result_type boost::math::tools::detail::continued_fraction_a_impl(Gen&, const U&, uintmax_t&) [with Gen = boost::math::detail::upper_incomplete_gamma_fract<double>; U = double; typename fraction_traits<Gen>::result_type = double; uintmax_t = long unsigned int]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/tools/fraction.hpp:311:44: required from ‘typename boost::math::tools::detail::fraction_traits<Gen>::result_type boost::math::tools::continued_fraction_a(Gen&, const U&) [with Gen = boost::math::detail::upper_incomplete_gamma_fract<double>; U = double; typename detail::fraction_traits<Gen>::result_type = double]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/gamma.hpp:378:68: required from ‘T boost::math::detail::upper_gamma_fraction(T, T, T) [with T = double]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/gamma.hpp:1679:44: required from ‘T boost::math::detail::gamma_incomplete_imp(T, T, bool, bool, const Policy&, T*) [with T = double; Policy = boost::math::policies::policy<boost::math::policies::pole_error<boost::math::policies::errno_on_error>, boost::math::policies::overflow_error<boost::math::policies::errno_on_error>, boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/gamma.hpp:2461:35: required from ‘boost::math::tools::promote_args_t<RT1, RT2> boost::math::gamma_p(RT1, RT2, const Policy&) [with RT1 = double; RT2 = double; Policy = policies::policy<policies::overflow_error<boost::math::policies::errno_on_error>, policies::pole_error<boost::math::policies::errno_on_error>, policies::promote_double<false>, policies::digits2<0>, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy>; tools::promote_args_t<RT1, RT2> = double]’
stan/lib/stan_math/stan/math/prim/fun/gamma_p.hpp:76:30: required from here
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/gamma.hpp:363:39: note: parameter passing for argument of type ‘std::pair<double, double>’ when C++17 is enabled changed to match C++14 in GCC 10.1
363 | BOOST_MATH_GPU_ENABLED result_type operator()()
| ^~~~~~~~In file included from stan/lib/stan_math/stan/math/prim/fun/owens_t.hpp:6,
from stan/lib/stan_math/stan/math/prim/fun.hpp:238:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/owens_t.hpp: In instantiation of ‘std::pair<_FIter, _FIter> boost::math::detail::owens_t_T1_accelerated(T, T, const Policy&) [with T = long double; Policy = boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/owens_t.hpp:874:46: required from ‘RealType boost::math::detail::owens_t_dispatch(RealType, RealType, RealType, const Policy&, const std::integral_constant<int, 65>&) [with RealType = long double; Policy = boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/owens_t.hpp:979:36: required from ‘RealType boost::math::detail::owens_t_dispatch(RealType, RealType, RealType, const Policy&) [with RealType = long double; Policy = boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/owens_t.hpp:1000:38: required from ‘RealType boost::math::detail::owens_t(RealType, RealType, const Policy&) [with RealType = long double; Policy = boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/owens_t.hpp:1072:86: required from ‘boost::math::tools::promote_args_t<RT1, RT2> boost::math::owens_t(T1, T2, const Policy&) [with T1 = double; T2 = double; Policy = policies::policy<policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy>; tools::promote_args_t<RT1, RT2> = double]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/owens_t.hpp:1078:24: required from ‘boost::math::tools::promote_args_t<RT1, RT2> boost::math::owens_t(T1, T2) [with T1 = double; T2 = double; tools::promote_args_t<RT1, RT2> = double]’
stan/lib/stan_math/stan/math/prim/fun/owens_t.hpp:62:30: required from here
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/owens_t.hpp:543:26: note: parameter passing for argument of type ‘std::pair<long double, long double>’ when C++17 is enabled changed to match C++14 in GCC 10.1
543 | std::pair<T, T> owens_t_T1_accelerated(T h, T a, const Policy& pol)
| ^~~~~~~~~~~~~~~~~~~~~~In file included from stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/binomial.hpp:17:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/beta.hpp: In instantiation of ‘boost::math::detail::ibeta_fraction2_t<T>::result_type boost::math::detail::ibeta_fraction2_t<T>::operator()() [with T = double; result_type = std::pair<double, double>]’:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/tools/fraction.hpp:134:20: required from ‘typename boost::math::tools::detail::fraction_traits<Gen>::result_type boost::math::tools::detail::continued_fraction_b_impl(Gen&, const U&, uintmax_t&) [with Gen = boost::math::detail::ibeta_fraction2_t<double>; U = double; typename fraction_traits<Gen>::result_type = double; uintmax_t = long unsigned int]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/tools/fraction.hpp:184:44: required from ‘typename boost::math::tools::detail::fraction_traits<Gen>::result_type boost::math::tools::continued_fraction_b(Gen&, const U&) [with Gen = boost::math::detail::ibeta_fraction2_t<double>; U = double; typename detail::fraction_traits<Gen>::result_type = double]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/beta.hpp:847:54: required from ‘T boost::math::detail::ibeta_fraction2(T, T, T, T, const Policy&, bool, T*) [with T = double; Policy = boost::math::policies::policy<boost::math::policies::pole_error<boost::math::policies::errno_on_error>, boost::math::policies::overflow_error<boost::math::policies::errno_on_error>, boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/beta.hpp:1499:36: required from ‘T boost::math::detail::ibeta_imp(T, T, T, const Policy&, bool, bool, T*) [with T = double; Policy = boost::math::policies::policy<boost::math::policies::pole_error<boost::math::policies::errno_on_error>, boost::math::policies::overflow_error<boost::math::policies::errno_on_error>, boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/beta.hpp:1536:20: required from ‘T boost::math::detail::ibeta_imp(T, T, T, const Policy&, bool, bool) [with T = double; Policy = boost::math::policies::policy<boost::math::policies::pole_error<boost::math::policies::errno_on_error>, boost::math::policies::overflow_error<boost::math::policies::errno_on_error>, boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/beta.hpp:1690:93: required from ‘boost::math::tools::promote_args_t<RT1, RT2, A> boost::math::ibeta(RT1, RT2, RT3, const Policy&) [with RT1 = double; RT2 = double; RT3 = double; Policy = policies::policy<policies::overflow_error<boost::math::policies::errno_on_error>, policies::pole_error<boost::math::policies::errno_on_error>, policies::promote_double<false>, policies::digits2<0>, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy>; tools::promote_args_t<RT1, RT2, A> = double]’
stan/lib/stan_math/stan/math/prim/fun/inc_beta.hpp:29:28: required from here
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/beta.hpp:810:39: note: parameter passing for argument of type ‘std::pair<double, double>’ when C++17 is enabled changed to match C++14 in GCC 10.1
810 | BOOST_MATH_GPU_ENABLED result_type operator()()
| ^~~~~~~~In file included from stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/gamma.hpp:19:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/tools/fraction.hpp: In function ‘typename boost::math::tools::detail::fraction_traits<Gen>::result_type boost::math::tools::detail::continued_fraction_a_impl(Gen&, const U&, uintmax_t&) [with Gen = boost::math::detail::upper_incomplete_gamma_fract<long double>; U = long double]’:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/tools/fraction.hpp:259:15: note: parameter passing for argument of type ‘std::pair<long double, long double>’ when C++17 is enabled changed to match C++14 in GCC 10.1
259 | value_type v = g();
| ^Start samplingRunning MCMC with 4 parallel chains...
Chain 1 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 1 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 1 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 1 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 1 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 1 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 1 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 1 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 1 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 1 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 1 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 2 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 2 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 2 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 2 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 2 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 2 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 2 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 2 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 2 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 2 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 3 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 3 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 3 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 3 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 3 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 3 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 3 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 3 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 3 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 3 Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4 Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4 Iteration: 100 / 2000 [ 5%] (Warmup)
Chain 4 Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4 Iteration: 300 / 2000 [ 15%] (Warmup)
Chain 4 Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4 Iteration: 500 / 2000 [ 25%] (Warmup)
Chain 4 Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4 Iteration: 700 / 2000 [ 35%] (Warmup)
Chain 4 Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4 Iteration: 900 / 2000 [ 45%] (Warmup)
Chain 4 Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4 Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4 Iteration: 1100 / 2000 [ 55%] (Sampling)
Chain 4 Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4 Iteration: 1300 / 2000 [ 65%] (Sampling)
Chain 4 Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4 Iteration: 1500 / 2000 [ 75%] (Sampling)
Chain 4 Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4 Iteration: 1700 / 2000 [ 85%] (Sampling)
Chain 4 Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4 Iteration: 1900 / 2000 [ 95%] (Sampling)
Chain 4 Iteration: 2000 / 2000 [100%] (Sampling)
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.1 seconds.Loading required namespace: rstan4.3.1 Summary
summary(normal_fit_brm) Family: gaussian
Links: mu = identity
Formula: y ~ 1
Data: normal_df (Number of observations: 100)
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 4.66 0.07 4.53 4.79 1.00 3560 2246
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.65 0.05 0.56 0.75 1.00 3307 2776
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).4.3.2 Draws (posterior samples)
normal_draws_brm <- posterior::as_draws_df(normal_fit_brm)
normal_draws_brm# A draws_df: 1000 iterations, 4 chains, and 5 variables
b_Intercept sigma Intercept lprior lp__
1 4.6 0.80 4.6 -4.4 -106
2 4.6 0.79 4.6 -4.4 -106
3 4.6 0.81 4.6 -4.4 -107
4 4.6 0.60 4.6 -4.4 -103
5 4.8 0.75 4.8 -4.4 -105
6 4.7 0.73 4.7 -4.4 -103
7 4.7 0.58 4.7 -4.4 -103
8 4.6 0.67 4.6 -4.4 -102
9 4.6 0.61 4.6 -4.4 -102
10 4.7 0.70 4.7 -4.4 -103
# ... with 3990 more draws
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}mcmc_trace(normal_draws_brm, pars = c("b_Intercept", "sigma"))
mcmc_hist(normal_draws_brm, pars = c("b_Intercept", "sigma"))`stat_bin()` using `bins = 30`. Pick better value `binwidth`.
4.3.3 Diagnostic summary
rstan::check_hmc_diagnostics(normal_fit_brm$fit)
Divergences:0 of 4000 iterations ended with a divergence.
Tree depth:0 of 4000 iterations saturated the maximum tree depth of 10.
Energy:E-BFMI indicated no pathological behavior.5 Poisson model
Count data (e.g., species richness, seed counts, pollinator visits…) often follow a Poisson distribution.
\[ y_i \sim \text{Poisson}(\lambda) \]
\[ y_i \sim \operatorname{Poisson\_log}(\log \lambda) \]
Q: What are the mean and variance (\(\lambda\)) of species richness per plot?
5.1 Dummy data
set.seed(123)
y <- rpois(n, lambda = 12.3)
pois_list <- list(
y = y,
N = n
)
pois_df <- tibble(
y = y,
N = n
)5.2 CmdStan
stan/pois.stan
data {
int<lower=0> N;
array[N] int<lower=0> y;
}
parameters {
real<lower=0> lambda;
}
model {
y ~ poisson(lambda);
lambda ~ normal(0, 10);
} stan/pois_re.stan
data {
int<lower=0> N;
array[N] int<lower=0> y;
}
parameters {
real log_lambda;
}
model {
y ~ poisson_log(log_lambda);
log_lambda ~ normal(0, 2.5);
} pois_mod <- cmdstan_model("stan/pois.stan")
pois_re_mod <- cmdstan_model("stan/pois_re.stan")pois_fit <- pois_mod$sample(
data = pois_list,
seed = 123,
chains = 4,
parallel_chains = 4,
iter_warmup = 1000, # number of warmup iterations
iter_sampling = 1000, # number of sampling iterations
refresh = 0 # print update every 500 iters
)Running MCMC with 4 parallel 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.2 seconds.pois_re_fit <- pois_re_mod$sample(
data = pois_list,
seed = 123,
chains = 4,
parallel_chains = 4,
iter_warmup = 1000, # number of warmup iterations
iter_sampling = 1000, # number of sampling iterations
refresh = 0 # print update every 500 iters
)Running MCMC with 4 parallel 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.2 seconds.pois_fit$summary()# A tibble: 2 × 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ 1836. 1836. 0.710 0.305 1834. 1836. 1.00 1898. 2556.
2 lambda 12.2 12.2 0.349 0.344 11.6 12.8 1.00 1370. 1871.pois_re_fit$summary()# A tibble: 2 × 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ 1833. 1834. 0.714 0.324 1.83e3 1.83e3 1.00 2170. 2574.
2 log_lambda 2.50 2.50 0.0287 0.0295 2.46e0 2.55e0 1.00 1367. 2197.5.3 brms
pois_fit_brm <- brm(y ~ 1,
family = poisson(),
data = pois_df,
refresh = 0,
backend = "cmdstanr")In file included from stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/hypergeometric_pFq.hpp:11,
from stan/lib/stan_math/stan/math/prim/fun/hypergeometric_pFq.hpp:8,
from stan/lib/stan_math/stan/math/prim/fun/hypergeometric_2F1.hpp:19,
from stan/lib/stan_math/stan/math/prim/fun/grad_2F1.hpp:14,
from stan/lib/stan_math/stan/math/prim/fun.hpp:117,
from stan/lib/stan_math/stan/math/prim/constraint.hpp:4,
from stan/lib/stan_math/stan/math/rev/constraint.hpp:24,
from stan/lib/stan_math/stan/math/rev.hpp:11,
from stan/lib/stan_math/stan/math.hpp:19,
from stan/src/stan/model/model_header.hpp:6:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp: In instantiation of ‘std::pair<_ForwardIterator, _ForwardIterator> boost::math::detail::hypergeometric_pFq_checked_series_impl(const Seq&, const Seq&, const Real&, const Policy&, const Terminal&, long long int&) [with Seq = std::vector<double>; Real = long double; Policy = boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>; Terminal = iteration_terminator]’:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/hypergeometric_pFq.hpp:61:107: required from ‘typename boost::math::tools::promote_args<Real, typename Seq::value_type>::type boost::math::hypergeometric_pFq(const Seq&, const Seq&, const Real&, Real*, const Policy&) [with Seq = std::vector<double>; Real = double; Policy = policies::policy<policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy>; typename tools::promote_args<Real, typename Seq::value_type>::type = double; typename Seq::value_type = double]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/hypergeometric_pFq.hpp:72:35: required from ‘typename boost::math::tools::promote_args<Real, typename Seq::value_type>::type boost::math::hypergeometric_pFq(const Seq&, const Seq&, const Real&, Real*) [with Seq = std::vector<double>; Real = double; typename tools::promote_args<Real, typename Seq::value_type>::type = double; typename Seq::value_type = double]’
stan/lib/stan_math/stan/math/prim/fun/hypergeometric_pFq.hpp:53:41: required from ‘stan::return_type_t<T, L, U> stan::math::hypergeometric_pFq(const Ta&, const Tb&, const Tz&) [with Ta = Eigen::Matrix<double, -1, 1>; Tb = Eigen::Matrix<double, -1, 1>; Tz = double; stan::require_all_vector_st<std::is_arithmetic, Ta, Tb>* <anonymous> = 0; stan::require_arithmetic_t<Tz>* <anonymous> = 0; stan::return_type_t<T, L, U> = double]’
stan/lib/stan_math/stan/math/prim/fun/hypergeometric_3F2.hpp:130:28: required from ‘auto stan::math::hypergeometric_3F2(const Ta&, const Tb&, const Tz&) [with Ta = std::vector<double>; Tb = std::vector<double>; Tz = double; stan::require_all_vector_t<T_y, T_high>* <anonymous> = 0; stan::require_stan_scalar_t<T2>* <anonymous> = 0]’
stan/lib/stan_math/stan/math/prim/fun/hypergeometric_3F2.hpp:154:28: required from ‘auto stan::math::hypergeometric_3F2(const std::initializer_list<_Tp>&, const std::initializer_list<_Value>&, const Tz&) [with Ta = double; Tb = double; Tz = double; stan::require_all_stan_scalar_t<T, L, U>* <anonymous> = 0]’
stan/lib/stan_math/stan/math/rev/fun/inv_inc_beta.hpp:80:37: required from here
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp:122:28: note: parameter passing for argument of type ‘std::pair<long double, long double>’ when C++17 is enabled changed to match C++14 in GCC 10.1
122 | std::pair<Real, Real> hypergeometric_pFq_checked_series_impl(const Seq& aj, const Seq& bj, const Real& z, const Policy& pol, const Terminal& termination, long long& log_scale)
| ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~In file included from stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/factorials.hpp:17,
from stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/binomial.hpp:16,
from stan/lib/stan_math/stan/math/prim/fun/choose.hpp:7,
from stan/lib/stan_math/stan/math/prim/functor/hcubature.hpp:21,
from stan/lib/stan_math/stan/math/prim/functor.hpp:14,
from stan/lib/stan_math/stan/math/rev/fun/to_arena.hpp:7,
from stan/lib/stan_math/stan/math/rev/core/operator_division.hpp:15,
from stan/lib/stan_math/stan/math/rev/core/operator_divide_equal.hpp:5,
from stan/lib/stan_math/stan/math/rev/core.hpp:29,
from stan/src/stan/model/model_base.hpp:8,
from stan/src/stan/model/model_header.hpp:4:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/gamma.hpp: In instantiation of ‘boost::math::detail::upper_incomplete_gamma_fract<T>::result_type boost::math::detail::upper_incomplete_gamma_fract<T>::operator()() [with T = double; result_type = std::pair<double, double>]’:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/tools/fraction.hpp:259:20: required from ‘typename boost::math::tools::detail::fraction_traits<Gen>::result_type boost::math::tools::detail::continued_fraction_a_impl(Gen&, const U&, uintmax_t&) [with Gen = boost::math::detail::upper_incomplete_gamma_fract<double>; U = double; typename fraction_traits<Gen>::result_type = double; uintmax_t = long unsigned int]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/tools/fraction.hpp:311:44: required from ‘typename boost::math::tools::detail::fraction_traits<Gen>::result_type boost::math::tools::continued_fraction_a(Gen&, const U&) [with Gen = boost::math::detail::upper_incomplete_gamma_fract<double>; U = double; typename detail::fraction_traits<Gen>::result_type = double]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/gamma.hpp:378:68: required from ‘T boost::math::detail::upper_gamma_fraction(T, T, T) [with T = double]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/gamma.hpp:1679:44: required from ‘T boost::math::detail::gamma_incomplete_imp(T, T, bool, bool, const Policy&, T*) [with T = double; Policy = boost::math::policies::policy<boost::math::policies::pole_error<boost::math::policies::errno_on_error>, boost::math::policies::overflow_error<boost::math::policies::errno_on_error>, boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/gamma.hpp:2461:35: required from ‘boost::math::tools::promote_args_t<RT1, RT2> boost::math::gamma_p(RT1, RT2, const Policy&) [with RT1 = double; RT2 = double; Policy = policies::policy<policies::overflow_error<boost::math::policies::errno_on_error>, policies::pole_error<boost::math::policies::errno_on_error>, policies::promote_double<false>, policies::digits2<0>, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy>; tools::promote_args_t<RT1, RT2> = double]’
stan/lib/stan_math/stan/math/prim/fun/gamma_p.hpp:76:30: required from here
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/gamma.hpp:363:39: note: parameter passing for argument of type ‘std::pair<double, double>’ when C++17 is enabled changed to match C++14 in GCC 10.1
363 | BOOST_MATH_GPU_ENABLED result_type operator()()
| ^~~~~~~~In file included from stan/lib/stan_math/stan/math/prim/fun/owens_t.hpp:6,
from stan/lib/stan_math/stan/math/prim/fun.hpp:238:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/owens_t.hpp: In instantiation of ‘std::pair<_FIter, _FIter> boost::math::detail::owens_t_T1_accelerated(T, T, const Policy&) [with T = long double; Policy = boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/owens_t.hpp:874:46: required from ‘RealType boost::math::detail::owens_t_dispatch(RealType, RealType, RealType, const Policy&, const std::integral_constant<int, 65>&) [with RealType = long double; Policy = boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/owens_t.hpp:979:36: required from ‘RealType boost::math::detail::owens_t_dispatch(RealType, RealType, RealType, const Policy&) [with RealType = long double; Policy = boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/owens_t.hpp:1000:38: required from ‘RealType boost::math::detail::owens_t(RealType, RealType, const Policy&) [with RealType = long double; Policy = boost::math::policies::policy<boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/owens_t.hpp:1072:86: required from ‘boost::math::tools::promote_args_t<RT1, RT2> boost::math::owens_t(T1, T2, const Policy&) [with T1 = double; T2 = double; Policy = policies::policy<policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy>; tools::promote_args_t<RT1, RT2> = double]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/owens_t.hpp:1078:24: required from ‘boost::math::tools::promote_args_t<RT1, RT2> boost::math::owens_t(T1, T2) [with T1 = double; T2 = double; tools::promote_args_t<RT1, RT2> = double]’
stan/lib/stan_math/stan/math/prim/fun/owens_t.hpp:62:30: required from here
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/owens_t.hpp:543:26: note: parameter passing for argument of type ‘std::pair<long double, long double>’ when C++17 is enabled changed to match C++14 in GCC 10.1
543 | std::pair<T, T> owens_t_T1_accelerated(T h, T a, const Policy& pol)
| ^~~~~~~~~~~~~~~~~~~~~~In file included from stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/binomial.hpp:17:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/beta.hpp: In instantiation of ‘boost::math::detail::ibeta_fraction2_t<T>::result_type boost::math::detail::ibeta_fraction2_t<T>::operator()() [with T = double; result_type = std::pair<double, double>]’:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/tools/fraction.hpp:134:20: required from ‘typename boost::math::tools::detail::fraction_traits<Gen>::result_type boost::math::tools::detail::continued_fraction_b_impl(Gen&, const U&, uintmax_t&) [with Gen = boost::math::detail::ibeta_fraction2_t<double>; U = double; typename fraction_traits<Gen>::result_type = double; uintmax_t = long unsigned int]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/tools/fraction.hpp:184:44: required from ‘typename boost::math::tools::detail::fraction_traits<Gen>::result_type boost::math::tools::continued_fraction_b(Gen&, const U&) [with Gen = boost::math::detail::ibeta_fraction2_t<double>; U = double; typename detail::fraction_traits<Gen>::result_type = double]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/beta.hpp:847:54: required from ‘T boost::math::detail::ibeta_fraction2(T, T, T, T, const Policy&, bool, T*) [with T = double; Policy = boost::math::policies::policy<boost::math::policies::pole_error<boost::math::policies::errno_on_error>, boost::math::policies::overflow_error<boost::math::policies::errno_on_error>, boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/beta.hpp:1499:36: required from ‘T boost::math::detail::ibeta_imp(T, T, T, const Policy&, bool, bool, T*) [with T = double; Policy = boost::math::policies::policy<boost::math::policies::pole_error<boost::math::policies::errno_on_error>, boost::math::policies::overflow_error<boost::math::policies::errno_on_error>, boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/beta.hpp:1536:20: required from ‘T boost::math::detail::ibeta_imp(T, T, T, const Policy&, bool, bool) [with T = double; Policy = boost::math::policies::policy<boost::math::policies::pole_error<boost::math::policies::errno_on_error>, boost::math::policies::overflow_error<boost::math::policies::errno_on_error>, boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy, boost::math::policies::default_policy>]’
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/beta.hpp:1690:93: required from ‘boost::math::tools::promote_args_t<RT1, RT2, A> boost::math::ibeta(RT1, RT2, RT3, const Policy&) [with RT1 = double; RT2 = double; RT3 = double; Policy = policies::policy<policies::overflow_error<boost::math::policies::errno_on_error>, policies::pole_error<boost::math::policies::errno_on_error>, policies::promote_double<false>, policies::digits2<0>, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy, policies::default_policy>; tools::promote_args_t<RT1, RT2, A> = double]’
stan/lib/stan_math/stan/math/prim/fun/inc_beta.hpp:29:28: required from here
stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/beta.hpp:810:39: note: parameter passing for argument of type ‘std::pair<double, double>’ when C++17 is enabled changed to match C++14 in GCC 10.1
810 | BOOST_MATH_GPU_ENABLED result_type operator()()
| ^~~~~~~~In file included from stan/lib/stan_math/lib/boost_1.87.0/boost/math/special_functions/gamma.hpp:19:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/tools/fraction.hpp: In function ‘typename boost::math::tools::detail::fraction_traits<Gen>::result_type boost::math::tools::detail::continued_fraction_a_impl(Gen&, const U&, uintmax_t&) [with Gen = boost::math::detail::upper_incomplete_gamma_fract<long double>; U = long double]’:
stan/lib/stan_math/lib/boost_1.87.0/boost/math/tools/fraction.hpp:259:15: note: parameter passing for argument of type ‘std::pair<long double, long double>’ when C++17 is enabled changed to match C++14 in GCC 10.1
259 | value_type v = g();
| ^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.5 seconds.summary(pois_fit_brm) Family: poisson
Links: mu = log
Formula: y ~ 1
Data: pois_df (Number of observations: 100)
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 2.50 0.03 2.44 2.56 1.00 1650 1850
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).6 Linear model
The simplest multiple linear regression model is the following, with two predictors (\(x_1\) and \(x_2\)), two slopes (\(\beta_1\) and \(\beta_2\)) and intercept coefficient (\(\beta_0\)), and normally distributed noise (\(\sigma\)). This model can be written using standard regression notation as
\[ y_i \sim N(\mu_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i}, \sigma) \]
or
\[ y_i \sim N(\mu_i = \boldsymbol{\beta \cdot x_i}, \sigma) \]
where \(\boldsymbol{\beta}\) is a 1 \(\times\) 3 coefficient matrix (or row vector) and \(\boldsymbol{x_i}\) is a 3 \(\times\) 1 predictor matrix (including the intercept). We use \(\boldsymbol{x}\) (3 \(\times\) N) to compute the \(\mu\) vector with length \(N\).
Q: What are the posterior means and 95% credible intervals of the coefficients (\(\beta_0\), \(\beta_1\), \(\beta_2\))?
set.seed(123)
n <- 100
beta <- c(2, 1.2, -0.8)
x1 <- rnorm(n, mean = 0, sd = 1)
x2 <- rnorm(n, mean = 0, sd = 1)
sigma <- 0.4
y <- rnorm(n, mean = beta[1] + beta[2] * x1 + beta[3] * x2, sd = sigma)lm_list <- list(
y = y,
N = length(y),
x = rbind(1, x1, x2) # 3 x N matrix
)
lm_df <- tibble(
y = y,
N = length(y),
x1 = x1,
x2 = x2
)6.1 CmdStan
stan/lm.stan
data {
int<lower=0> N;
vector[N] y;
matrix[3, N] x;
}
parameters {
row_vector[3] beta;
real<lower=0> sigma;
}
model {
y ~ normal(beta * x, sigma);
beta ~ normal(0, 2.5);
sigma ~ cauchy(0, 2.5);
} lm_mod <- cmdstan_model("stan/lm.stan")
lm_fit <- lm_mod$sample(
data = lm_list,
seed = 123,
chains = 4,
parallel_chains = 4,
iter_warmup = 1000, # number of warmup iterations
iter_sampling = 1000, # number of sampling iterations
refresh = 0 # don't print update
)Running MCMC with 4 parallel 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.2 seconds.Use posterior::summarise_draws to check more detailed summary statistics.
my_summary <- function(draws, probs = c(0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975)) {
posterior::summarise_draws(
draws,
mean = ~mean(.x),
sd = ~sd(.x),
~posterior::quantile2(.x, probs = probs)
)
}
my_summary(lm_fit$draws())# A tibble: 5 × 10
variable mean sd q2.5 q5 q25 q50 q75 q95 q97.5
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ 44.6 1.43 41.0 41.9 43.9 44.9 45.6 46.2 46.4
2 beta[1] 2.05 0.0393 1.98 1.99 2.03 2.05 2.08 2.12 2.13
3 beta[2] 1.15 0.0417 1.06 1.08 1.12 1.15 1.17 1.21 1.23
4 beta[3] -0.791 0.0410 -0.869 -0.857 -0.819 -0.791 -0.763 -0.724 -0.708
5 sigma 0.385 0.0285 0.333 0.341 0.365 0.383 0.403 0.433 0.4446.2 brms
lm_fit_brm <- brm(
y ~ x1 + x2,
family = gaussian(),
data = lm_df,
prior = c(
prior(normal(0, 2.5), class = "Intercept"),
prior(normal(0, 2.5), class = "b"),
prior(cauchy(0, 2.5), class = "sigma")
),
refresh = 0, # don't print update
backend = "cmdstanr"
)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.0 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.0 seconds.
Total execution time: 0.5 seconds.summary(lm_fit_brm) Family: gaussian
Links: mu = identity
Formula: y ~ x1 + x2
Data: lm_df (Number of observations: 100)
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 2.05 0.04 1.98 2.13 1.00 4619 3168
x1 1.15 0.04 1.06 1.23 1.00 4356 2728
x2 -0.79 0.04 -0.87 -0.71 1.00 4392 3066
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.39 0.03 0.33 0.45 1.00 4206 2649
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).7 Varying intercepts
H: There is a power-law relationship (\(y =\beta_0x^{\beta_1}\)) between tree diameter (DBH) and tree maximum height, and the scaling factor \(\beta_0\) varies among species.
\[ \text{log}\; y_{ij} \sim N(\text{log}\; \beta_{0j} + \beta_1 x_{i}, \sigma) \]
\[ \text{log}\; \beta_{0j} \sim N(\mu_0, \tau) \]
\[ \mu_0 \sim N(0, 2.5) \]
\[ \tau \sim \text{Half-Cauchy}(0, 1) \]
Q: What are the species-level scaling factors (\(\beta_{0j}\)) and the global scaling exponent (\(\beta_1\))?
7.1 Dummy data
- We simulate a dataset with 10 species and 20 trees per species.
- Each species has wood density (wd) values.
set.seed(12345)
n_sp <- 10
n_rep <- 20
wd <- rnorm(n_sp, 0, 1)
gamma0 <- 0.6
gamma1 <- 0.1
sigma_y <- 0.1
b1_hat <- gamma1 * wd + gamma0
b1 <- rnorm(n_sp, b1_hat, 0.01)
log_b0 <- rnorm(n_sp, 0.55, 0.05)
# ---- simulate ----
allo_df <- tibble(
sp = factor(rep(paste0("sp", 1:n_sp), each = n_rep)),
wd = rep(wd, each = n_rep),
# now log_xx ~ Normal(mean log-dbh, sd = 0.5)
log_xx = rnorm(n_sp * n_rep, mean = log(40), sd = 0.5)) |>
mutate(
# add observation-level noise on log-height
log_y = rnorm(
n(),
log_b0[as.integer(sp)] + b1[as.integer(sp)] * log_xx,
sigma_y),
dbh = exp(log_xx),
h = exp(log_y)) |>
select(sp, wd, dbh, h)
dbh_hist <- allo_df |>
ggplot(aes(dbh)) +
geom_histogram() +
xlab("DBH (cm)") +
theme(
legend.position = "none",
plot.title = element_text(size = 24)
)
h_hist <- allo_df |>
ggplot(aes(h)) +
geom_histogram() +
xlab("Height (m)") +
theme(
legend.position = "none",
plot.title = element_text(size = 24)
)
p1 <- allo_df |>
ggplot(aes(x = dbh, y = h, col = sp)) +
geom_point() +
xlab("DBH (cm)") +
ylab(expression("Height (m)")) +
theme(
legend.position = "none",
plot.title = element_text(size = 24)
)
p2 <- p1 +
scale_x_log10() +
scale_y_log10()
dbh_hist + h_hist + p1 + p2
allo_list <- list(
log_h = log(allo_df$h),
log_dbh = log(allo_df$dbh),
sp = as.integer(allo_df$sp),
N = nrow(allo_df),
J = allo_df$sp |> unique() |> length()
)7.2 CmdStan
stan/vint.stan
data {
int<lower=1> N; // total number of observations
int<lower=1> J; // number of species
array[N] int<lower=1, upper=J> sp; // sp index
vector[N] log_dbh; // predictor: log(DBH)
vector[N] log_h; // response: log(height)
}
parameters {
real mu0; // grand mean of log(mu₀)
real<lower=0> tau; // SD of species intercepts
vector[J] log_beta0; // species‐level intercepts = log(β₀ⱼ)
real beta1; // common slope on log(DBH)
real<lower=0> sigma; // residual SD
}
model {
// Likelihood
log_h ~ normal(log_beta0[sp] + beta1 * log_dbh, sigma);
// Priors
mu0 ~ normal(0, 2.5);
tau ~ cauchy(0, 1);
log_beta0 ~ normal(mu0, tau);
beta1 ~ normal(0, 2.5);
sigma ~ cauchy(0, 1);
}
// we need this for model selection
generated quantities {
vector[N] log_lik;
for (i in 1:N)
log_lik[i] = normal_lpdf(log_h[i]
| log_beta0[sp[i]] + beta1 * log_dbh[i],
sigma);
} vint_mod <- cmdstan_model("stan/vint.stan")
vint_fit <- vint_mod$sample(
data = allo_list,
seed = 123,
chains = 4,
parallel_chains = 4,
iter_warmup = 1000, # number of warmup iterations
iter_sampling = 1000, # number of sampling iterations
refresh = 0 # don't print update
)Running MCMC with 4 parallel chains...Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp10tuT6/model-23ced1a250.stan', line 19, column 2 to column 57)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 2 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp10tuT6/model-23ced1a250.stan', line 19, column 2 to column 57)Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 2 Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 3 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmp10tuT6/model-23ced1a250.stan', line 19, column 2 to column 57)Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 3 Chain 1 finished in 0.3 seconds.
Chain 2 finished in 0.3 seconds.
Chain 3 finished in 0.3 seconds.
Chain 4 finished in 0.3 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.3 seconds.
Total execution time: 0.4 seconds.#vint_fit$summary()
my_summary(vint_fit$draws())# A tibble: 215 × 10
variable mean sd q2.5 q5 q25 q50 q75 q95
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ 365. 2.71 359. 360. 364. 366. 3.67e+2 3.69e+2
2 mu0 0.479 0.121 0.239 0.280 0.403 0.479 5.54e-1 6.68e-1
3 tau 0.333 0.0909 0.208 0.220 0.269 0.315 3.79e-1 5.05e-1
4 log_beta0[… 0.779 0.0565 0.671 0.686 0.740 0.778 8.16e-1 8.74e-1
5 log_beta0[… 0.918 0.0542 0.812 0.830 0.880 0.916 9.55e-1 1.01e+0
6 log_beta0[… 0.436 0.0554 0.331 0.345 0.398 0.436 4.74e-1 5.28e-1
7 log_beta0[… 0.261 0.0582 0.150 0.166 0.222 0.261 3.01e-1 3.58e-1
8 log_beta0[… 0.617 0.0534 0.514 0.529 0.580 0.617 6.53e-1 7.06e-1
9 log_beta0[… -0.0183 0.0546 -0.124 -0.107 -0.0563 -0.0185 1.93e-2 7.34e-2
10 log_beta0[… 0.707 0.0534 0.605 0.622 0.669 0.706 7.43e-1 7.96e-1
# ℹ 205 more rows
# ℹ 1 more variable: q97.5 <dbl>vint_fit$diagnostic_summary()$num_divergent
[1] 0 0 0 0
$num_max_treedepth
[1] 0 0 0 0
$ebfmi
[1] 1.0810518 0.9796202 1.0400008 0.99934187.3 brms
priors <- c(
prior(normal(0, 2.5), class = Intercept),
prior(normal(0, 2.5), class = b),
prior(cauchy(0, 1), class = sd),
prior(cauchy(0, 1), class = sigma)
)
vint_fit_brm <- brm(
log(h) ~ log(dbh) + (1 | sp),
family = gaussian(),
data = allo_df,
prior = priors,
seed = 123,
refresh = 0, # don't print update
backend = "cmdstanr"
)Running MCMC with 4 sequential chains...
Chain 1 finished in 0.7 seconds.
Chain 2 finished in 0.7 seconds.
Chain 3 finished in 0.7 seconds.
Chain 4 finished in 0.6 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.7 seconds.
Total execution time: 2.9 seconds.summary(vint_fit_brm) Family: gaussian
Links: mu = identity
Formula: log(h) ~ log(dbh) + (1 | sp)
Data: allo_df (Number of observations: 200)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Multilevel Hyperparameters:
~sp (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.34 0.10 0.21 0.58 1.00 590 1242
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.48 0.12 0.25 0.71 1.00 897 1238
logdbh 0.61 0.01 0.58 0.64 1.00 2371 2261
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.10 0.01 0.09 0.11 1.00 2165 2143
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).plot(vint_fit_brm)
7.4 lme4
vint_fit_lme <- lme4::lmer(
log(h) ~ log(dbh) + (1 | sp),
data = allo_df)
summary(vint_fit_lme)Linear mixed model fit by REML ['lmerMod']
Formula: log(h) ~ log(dbh) + (1 | sp)
Data: allo_df
REML criterion at convergence: -298.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.66966 -0.62316 0.04254 0.60122 2.28315
Random effects:
Groups Name Variance Std.Dev.
sp (Intercept) 0.084434 0.29057
Residual 0.009794 0.09897
Number of obs: 200, groups: sp, 10
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.47916 0.10454 4.584
log(dbh) 0.61154 0.01314 46.525
Correlation of Fixed Effects:
(Intr)
log(dbh) -0.4727.5 Visualization and predictions
7.5.1 CmdStan
# Predicted DBH-height curves per species with 95% intervals
pred_grid <- allo_df |>
dplyr::group_by(sp) |>
dplyr::summarise(
dbh = list(seq(min(dbh), max(dbh), length.out = 60)),
.groups = "drop"
) |>
tidyr::unnest(dbh)
cmd_draws <- posterior::as_draws_matrix(vint_fit$draws(c("beta1", "log_beta0")))
beta1_draws <- cmd_draws[, "beta1"]
log_beta0_cols <- grep("^log_beta0\\[", colnames(cmd_draws), value = TRUE)
species_levels <- levels(allo_df$sp)DBH for predicted lines for each species.
pred_grid# A tibble: 600 × 2
sp dbh
<fct> <dbl>
1 sp1 12.2
2 sp1 14.0
3 sp1 15.8
4 sp1 17.6
5 sp1 19.5
6 sp1 21.3
7 sp1 23.1
8 sp1 25.0
9 sp1 26.8
10 sp1 28.6
# ℹ 590 more rowsPosterior draws of the slope parameter (\(\beta_1\)).
str(beta1_draws) 'draws_matrix' num [1:4000, 1] 0.615 0.614 0.598 0.609 0.61 ...
- attr(*, "dimnames")=List of 2
..$ draw : chr [1:4000] "1" "2" "3" "4" ...
..$ variable: chr "beta1"
- attr(*, "nchains")= int 4Column names for log-transformed species-level intercepts (\(\log \beta_{0j}\)). We use this to extract posterior draws for each species.
str(log_beta0_cols) chr [1:10] "log_beta0[1]" "log_beta0[2]" "log_beta0[3]" "log_beta0[4]" ...pred_cmd_summary <- dplyr::bind_cols(
pred_grid,
purrr::map2_dfr(
match(pred_grid$sp, species_levels),
log(pred_grid$dbh),
~{
log_mu <- cmd_draws[, log_beta0_cols[.x]] + beta1_draws * .y
height_draws <- exp(log_mu)
tibble::tibble(
estimate = stats::median(height_draws),
lower = stats::quantile(height_draws, 0.025),
upper = stats::quantile(height_draws, 0.975)
)
}
)
)Compute predicted height summary for each species and DBH.
pred_cmd_summary# A tibble: 600 × 5
sp dbh estimate lower upper
<fct> <dbl> <dbl> <dbl> <dbl>
1 sp1 12.2 10.0 9.50 10.7
2 sp1 14.0 10.9 10.4 11.6
3 sp1 15.8 11.8 11.2 12.5
4 sp1 17.6 12.6 12.0 13.3
5 sp1 19.5 13.4 12.8 14.1
6 sp1 21.3 14.2 13.5 14.9
7 sp1 23.1 14.9 14.2 15.6
8 sp1 25.0 15.6 14.9 16.4
9 sp1 26.8 16.3 15.6 17.1
10 sp1 28.6 17.0 16.2 17.8
# ℹ 590 more rowsggplot(pred_cmd_summary, aes(dbh, estimate, colour = sp)) +
geom_ribbon(aes(ymin = lower, ymax = upper, fill = sp), alpha = 0.15, colour = NA) +
geom_line() +
geom_point(
data = allo_df,
aes(dbh, h, colour = sp),
size = 1,
alpha = 0.6,
inherit.aes = FALSE
) +
scale_x_log10() +
scale_y_log10() +
guides(fill = "none") +
labs(
x = "DBH (cm)",
y = "Height (m)",
colour = "Species"
)
7.5.2 brms
With fitted() available, preparing pred_summary is simpler than in the CmdStan workflow.
pred_summary <- fitted(vint_fit_brm, newdata = pred_grid, re_formula = NULL) |>
tibble::as_tibble() |>
dplyr::bind_cols(pred_grid) |>
dplyr::mutate(
estimate = exp(Estimate),
lower = exp(Q2.5),
upper = exp(Q97.5)
)
ggplot(pred_summary, aes(dbh, estimate, color = sp)) +
geom_ribbon(aes(ymin = lower, ymax = upper, fill = sp), alpha = 0.15, colour = NA) +
geom_line() +
geom_point(
data = allo_df,
aes(dbh, h, colour = sp),
size = 1,
alpha = 0.6,
inherit.aes = FALSE
) +
scale_x_log10() +
scale_y_log10() +
guides(fill = "none") +
labs(
x = "DBH (cm)",
y = "Height (m)",
colour = "Species"
)
8 Reparameterization for multilevel models
Diagnosing Biased Inference with Divergences by Michael Betancourt
8.1 Eight schools example
What is the effect of the treatment (a new education method) on students’ scores for each school and across schools?
| school | y | sigma |
|---|---|---|
| School_1 | 28 | 15 |
| School_2 | 8 | 10 |
| School_3 | -3 | 16 |
| School_4 | 7 | 11 |
| School_5 | -1 | 9 |
| School_6 | 1 | 11 |
| School_7 | 18 | 10 |
| School_8 | 12 | 18 |
8.2 Centered parameterization
\[ y_j \sim \mathcal{N}\bigl(\theta_j,\;\sigma_j\bigr) \]
\[ \theta_j \sim \mathcal{N}\bigl(\mu,\;\tau \bigr) \]
\[ \mu \sim \mathcal{N}\bigl(0,\;5 \bigr) \]
\[ \tau \sim \text{Half-Cauchy}(0, 2.5) \]
Note: \(\sigma_j\) is known constant (i.e., data), we don’t have to estimate it.
This parameterization is intuitive, but it often leads to convergence issues in hierarchical models.
8.3 Non-centered parameterization
\[ \tilde{\theta_j} \sim \mathcal{N}\bigl(0,\;1 \bigr) \]
\[ \theta_j = \mu + \tau \cdot \tilde{\theta_j} \]
In this parameterization, we introduce a new latent variable \(\tilde{\theta_j}\) that is independent of the other parameters. This allows the sampler to explore the posterior more efficiently and avoids problems with convergence. We often use this parameterization when we have hierarchical models with varying intercepts or slopes.
9 Varying slopes and intercepts
H: There is a power-law relationship (\(y =\beta_0 x^{\beta_1}\)) between tree diameter (DBH) and tree maximum height, and both the scaling factor \(\beta_0\) and the exponent \(\beta_1\) vary among species.
\[ \text{log}\; y_{ij} \sim N(\text{log}\; \beta_{0j} + \beta_{1j} x_{ij}, \sigma) \]
\[ \begin{bmatrix} \text{log}\; \beta_{0j} \\ \beta_{1j} \end{bmatrix} \sim \mathrm{MVN} \Biggl( \boldsymbol{\mu} = \begin{bmatrix} \mu_{0} \\ \mu_{1} \end{bmatrix}, \boldsymbol{\Sigma} = \begin{bmatrix} \sigma_{0}^2 & \rho \sigma_{0} \sigma_{1} \\ \rho \sigma_{0} \sigma_{1} & \sigma_{1}^2 \end{bmatrix} \Biggr) \]
For reference, the multivariate normal (MVN) distribution looks like this:
Q: What are the species-level scaling factors (\(\beta_{0j}\)) and scaling exponent (\(\beta_{0j}\))? What are the global scaling exponent (\(\mu_0\)) and scaling exponent (\(\mu_1\))?
9.1 Centered parameterization
\[ \boldsymbol{\mu} \sim N(0, 2.5) \]
\[ \boldsymbol{\Sigma} \sim \mathcal{W}^{-1}_p(\nu,\,\Psi) \]
The inverse-whishart distribution \(\mathcal{W}^{-1}\) is a conjugate for the covariance matrix \(\boldsymbol{\Sigma}\) in MVN. However, despite conjugacy, inverse–Wishart priors are often computationally difficult to sample.
Recommendation: use a Cholesky factorization with independent priors on scales and correlation (below).
9.2 Non-Centered parameterization
Optimization through Cholesky factorization
\[ \boldsymbol{z}_j \sim \mathcal{N}(\mathbf0,\,I_2) \] \[ \begin{pmatrix}\log\beta_{0j}\\[3pt]\beta_{1j}\end{pmatrix} = \boldsymbol{\mu} + L\, \boldsymbol{z}_j \]
Some linear algebra (Cholesky decomposition):
\[ \Sigma \;=\;L\,L^\top \;=\;\bigl[\mathrm{diag}(\tau)\,L_\Omega\bigr]\, \bigl[L_\Omega^\top\,\mathrm{diag}(\tau)\bigr] \;=\;\mathrm{diag}(\tau)\,\Omega\,\mathrm{diag}(\tau). \]
\[ L_\Omega \sim \mathrm{LKJ\_Corr\_Cholesky}(2) \quad \text{or} \quad \Omega \sim \mathrm{LKJ\_Corr}(2) \] \[ \tau_i \sim \mathrm{Half\text{-}Cauchy}(0,1) \]
Rather than sampling the covariance matrix \(\Sigma\) or correlation matrix \(\Omega\) directly, sampling the Cholesky factor \(L_\Omega\) is computationally more efficient and numerically stable.
In the LKJ() prior:
\(\eta = 1\) indicates a flat prior (uniform distribution) over correlation matrices.
\(\eta > 1\) indicates weaker correlations.
We usually use \(\eta = 2\) because we assume at least weak correlations among the parameters while still allowing moderate to strong correlations.
9.3 CmdStan
stan/vslope.stan
data {
int<lower=0> N; // num individuals
int<lower=1> K; // num ind predictors
int<lower=1> J; // num groups
array[N] int<lower=1, upper=J> jj; // group for individual
int<lower=1> L; // number of group-level predictors
matrix[N, K] x; // individual predictors
matrix[L, J] u; // group predictors
vector[N] y; // outcomes
}
parameters {
matrix[K, L] gamma; // coefficients across groups
real<lower=0> sigma; // prediction error scale
vector<lower=0>[K] tau; // prior scale
matrix[K, J] z;
cholesky_factor_corr[K] L_Omega;
}
transformed parameters {
matrix[K, J] beta;
// mu = gamma * u
beta = gamma * u + diag_pre_multiply(tau, L_Omega) * z;
}
model {
to_vector(z) ~ std_normal();
tau ~ cauchy(0, 2.5);
L_Omega ~ lkj_corr_cholesky(2);
to_vector(gamma) ~ normal(0, 2.5);
for (n in 1:N) {
y[n] ~ normal(x[n, ] * beta[, jj[n]], sigma);
}
} allo_vslope_list <- list(
y = log(allo_df$h),
x = cbind(1, log(allo_df$dbh)),
jj = as.integer(allo_df$sp),
N = nrow(allo_df),
K = 2, # number of predictors (intercept + slope)
J = allo_df$sp |> unique() |> length(),
L = 1 # number of group-level predictors (intercept)
)
allo_vslope_list$u <- matrix(1, nrow = allo_vslope_list$L, ncol = allo_vslope_list$J)vslope_mod <- cmdstan_model("stan/vslope.stan")
vslope_fit <- vslope_mod$sample(
data = allo_vslope_list,
seed = 123,
chains = 4,
parallel_chains = 4,
iter_warmup = 1000, # number of warmup iterations
iter_sampling = 1000, # number of sampling iterations
adapt_delta = 0.99, # increase adapt_delta to avoid divergent transitions
max_treedepth = 15, # increase max_treedepth to avoid max treedepth errors
refresh = 0, # don't print update
show_messages = FALSE # suppress Stan’s “Informational Message” output
)Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 32, column 4 to column 49)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 2 Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 2 Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 2 Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 2 Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 2 Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 3 Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 3 Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 3 Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 3 Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 3 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 my_summary(vslope_fit$draws())# A tibble: 50 × 10
variable mean sd q2.5 q5 q25 q50 q75 q95
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ 3.60e+2 4.39 3.50e+2 3.52e+2 357. 3.60e+2 3.63e+2 366.
2 gamma[1,1] 5.01e-1 0.0555 3.91e-1 4.10e-1 0.465 5.02e-1 5.38e-1 0.591
3 gamma[2,1] 6.07e-1 0.0313 5.47e-1 5.56e-1 0.587 6.06e-1 6.26e-1 0.660
4 sigma 9.24e-2 0.00481 8.37e-2 8.49e-2 0.0889 9.21e-2 9.55e-2 0.101
5 tau[1] 7.10e-2 0.0570 2.34e-3 4.79e-3 0.0278 5.94e-2 9.96e-2 0.181
6 tau[2] 8.75e-2 0.0290 4.73e-2 5.14e-2 0.0679 8.24e-2 1.01e-1 0.142
7 z[1,1] 6.05e-1 0.895 -1.27e+0 -9.15e-1 0.0410 6.46e-1 1.23e+0 2.01
8 z[2,1] 6.61e-1 0.653 -7.88e-1 -5.21e-1 0.286 7.03e-1 1.09e+0 1.65
9 z[1,2] 1.95e-1 0.975 -1.83e+0 -1.50e+0 -0.441 2.23e-1 8.80e-1 1.73
10 z[2,2] 1.30e+0 0.658 -4.72e-2 1.93e-1 0.900 1.32e+0 1.73e+0 2.33
# ℹ 40 more rows
# ℹ 1 more variable: q97.5 <dbl>vslope_fit$diagnostic_summary()$num_divergent
[1] 0 0 0 0
$num_max_treedepth
[1] 0 0 0 0
$ebfmi
[1] 0.8136123 0.8535260 0.8342642 0.76026239.4 brms
priors <- c(
prior(normal(0, 2.5), class = Intercept),
prior(normal(0, 2.5), class = b),
prior(lkj(2), class = cor), # Ω ~ LKJ(2)
prior(cauchy(0, 2.5), class = sigma)
)
slope_fit_brm <- brm(
log(h) ~ log(dbh) + (1 + log(dbh) | sp),
family = gaussian(),
data = allo_df,
prior = priors,
refresh = 0, # don't print update
seed = 123,
chains = 4, cores = 4,
control = list(
adapt_delta = 0.99, # default is 0.8
max_treedepth = 15 # default is 10
),
backend = "cmdstanr"
)Running MCMC with 4 parallel chains...
Chain 1 finished in 7.2 seconds.
Chain 4 finished in 9.3 seconds.
Chain 3 finished in 10.0 seconds.
Chain 2 finished in 10.3 seconds.
All 4 chains finished successfully.
Mean chain execution time: 9.2 seconds.
Total execution time: 10.4 seconds.summary(slope_fit_brm) Family: gaussian
Links: mu = identity
Formula: log(h) ~ log(dbh) + (1 + log(dbh) | sp)
Data: allo_df (Number of observations: 200)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Multilevel Hyperparameters:
~sp (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(Intercept) 0.07 0.06 0.00 0.21 1.00 1853
sd(logdbh) 0.09 0.03 0.05 0.16 1.00 1320
cor(Intercept,logdbh) 0.09 0.42 -0.74 0.84 1.00 672
Tail_ESS
sd(Intercept) 1523
sd(logdbh) 1835
cor(Intercept,logdbh) 1251
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.50 0.05 0.40 0.61 1.00 3410 3243
logdbh 0.60 0.03 0.54 0.66 1.00 1362 1765
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.09 0.00 0.08 0.10 1.00 6219 2878
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).9.5 lme4
slope_fit_lme <- lme4::lmer(
log(h) ~ log(dbh) + (1 + log(dbh) | sp),
data = allo_df)boundary (singular) fit: see help('isSingular')summary(slope_fit_lme)Linear mixed model fit by REML ['lmerMod']
Formula: log(h) ~ log(dbh) + (1 + log(dbh) | sp)
Data: allo_df
REML criterion at convergence: -327.2
Scaled residuals:
Min 1Q Median 3Q Max
-3.14307 -0.63987 0.08714 0.61696 3.03735
Random effects:
Groups Name Variance Std.Dev. Corr
sp (Intercept) 0.001326 0.03641
log(dbh) 0.004607 0.06787 1.00
Residual 0.008397 0.09164
Number of obs: 200, groups: sp, 10
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.49980 0.04759 10.50
log(dbh) 0.60637 0.02468 24.57
Correlation of Fixed Effects:
(Intr)
log(dbh) -0.264
optimizer (nloptwrap) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')9.6 Visualization
Almost the same as before, but now we have varying slopes as well as intercepts.
9.6.1 CmdStan
vslope_matrix <- posterior::as_draws_matrix(vslope_fit$draws("beta"))
intercept_cols <- grep("^beta\\[1,", colnames(vslope_matrix), value = TRUE)
slope_cols <- grep("^beta\\[2,", colnames(vslope_matrix), value = TRUE)
str(vslope_matrix) 'draws_matrix' num [1:4000, 1:20] 0.749 0.482 0.469 0.633 0.688 ...
- attr(*, "dimnames")=List of 2
..$ draw : chr [1:4000] "1" "2" "3" "4" ...
..$ variable: chr [1:20] "beta[1,1]" "beta[2,1]" "beta[1,2]" "beta[2,2]" ...
- attr(*, "nchains")= int 4str(intercept_cols) chr [1:10] "beta[1,1]" "beta[1,2]" "beta[1,3]" "beta[1,4]" "beta[1,5]" ...str(slope_cols) chr [1:10] "beta[2,1]" "beta[2,2]" "beta[2,3]" "beta[2,4]" "beta[2,5]" ...pred_grid_vslope <- pred_grid |>
dplyr::mutate(log_dbh = log(dbh))
pred_cmd_summary <- dplyr::bind_cols(
pred_grid_vslope,
purrr::map2_dfr(
match(pred_grid_vslope$sp, species_levels),
pred_grid_vslope$log_dbh,
~{
intercept_draws <- vslope_matrix[, intercept_cols[.x]]
slope_draws <- vslope_matrix[, slope_cols[.x]]
log_mu <- intercept_draws + slope_draws * .y
height_draws <- exp(log_mu)
tibble::tibble(
estimate = stats::median(height_draws),
lower = stats::quantile(height_draws, 0.025),
upper = stats::quantile(height_draws, 0.975)
)
}
)
)
ggplot(pred_cmd_summary, aes(dbh, estimate, colour = sp)) +
geom_ribbon(
aes(ymin = lower, ymax = upper, fill = sp),
alpha = 0.15,
colour = NA
) +
geom_line() +
geom_point(
data = allo_df,
aes(dbh, h, colour = sp),
size = 1,
alpha = 0.6,
inherit.aes = FALSE
) +
scale_x_log10() +
scale_y_log10() +
guides(fill = "none", colour = guide_legend(title = "Species")) +
labs(
x = "Diameter at breast height (cm)",
y = "Height (m)",
title = "Posterior predictive curves from CmdStan model"
)
9.6.2 brms
pred_brm_summary <- fitted(
slope_fit_brm,
newdata = pred_grid,
re_formula = NULL
) |>
tibble::as_tibble() |>
dplyr::bind_cols(pred_grid) |>
dplyr::mutate(
estimate = exp(Estimate),
lower = exp(Q2.5),
upper = exp(Q97.5)
)
ggplot(pred_brm_summary, aes(dbh, estimate, colour = sp)) +
geom_ribbon(
aes(ymin = lower, ymax = upper, fill = sp),
alpha = 0.15,
colour = NA
) +
geom_line() +
geom_point(
data = allo_df,
aes(dbh, h, colour = sp),
size = 1,
alpha = 0.6,
inherit.aes = FALSE
) +
scale_x_log10() +
scale_y_log10() +
guides(fill = "none", colour = guide_legend(title = "Species")) +
labs(
x = "Diameter at breast height (cm)",
y = "Height (m)",
title = "Posterior predictive curves from brms model"
)
10 Varying slopes and intercepts, and group-level predictors
H: There is a power-law relationship (\(y =\beta_0 x^{\beta_1}\)) between tree diameter (DBH) and tree maximum height, and both the scaling factor \(\beta_0\) and the exponent \(\beta_1\) vary among species. Those species-level variations can be explained by wood density (group-level predictor: \(u\)) of each species.
\[ \text{log}\; y_{ij} \sim N(\text{log}\; \beta_{0j} + \beta_{1j} x_{ij}, \sigma) \]
\[ \begin{bmatrix} \text{log}\; \beta_{0j} \\ \beta_{1j} \end{bmatrix} \sim \mathrm{MVN} \Biggl( \underbrace{ \boldsymbol{\mu} = \begin{bmatrix} \gamma_{0}^{(\beta_0)} + \gamma_{1}^{(\beta_0)}u_{j} \\ \gamma_{0}^{(\beta_1)} + \gamma_{1}^{(\beta_1)}u_{j} \end{bmatrix}}_{\text{mean depends on }u_j}, \boldsymbol{\Sigma} = \begin{bmatrix} \sigma_{0}^2 & \rho \sigma_{0} \sigma_{1} \\ \rho \sigma_{0} \sigma_{1} & \sigma_{1}^2 \end{bmatrix} \Biggr) \]
Q: What are the coefficients for species-level regressions (\(\gamma_{0}^{(\beta_0)}\), \(\gamma_{1}^{(\beta_0}\), \(\gamma_{0}^{(\beta_1)}\), and \(\gamma_{1}^{(\beta_1}\))?
10.1 CmdStan
Now we have two group-level predictors (intercept and slope), and \(u\) includes the intercept and wood density. The rest of this section follows the setup from the varying slopes and intercepts section.
allo_vslope_list$L <- 2
allo_vslope_list$u <- rbind(
rep(1, allo_vslope_list$J),
trait
)
str(allo_vslope_list)List of 8
$ y : num [1:200] 3.29 3.66 3.7 3.58 3.19 ...
$ x : num [1:200, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
$ jj: int [1:200] 1 1 1 1 1 1 1 1 1 1 ...
$ N : int 200
$ K : num 2
$ J : int 10
$ L : num 2
$ u : num [1:2, 1:10] 1 0.586 1 0.709 1 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "" "trait"
.. ..$ : NULLvgrp_fit <- vslope_mod$sample(
data = allo_vslope_list,
seed = 1234,
chains = 4,
parallel_chains = 4,
iter_warmup = 1000, # number of warmup iterations
iter_sampling = 1000, # number of sampling iterations
adapt_delta = 0.95, # increase adapt_delta to avoid divergent transitions
max_treedepth = 15, # increase max_treedepth to avoid max treedepth errors
refresh = 0 # don't print update
)Running MCMC with 4 parallel chains...Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 1 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 1 Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 2 Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 2 Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 2 Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 2 Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 2 Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 2 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 2 Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 3 Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 3 Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 3 Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 3 Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 3 Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 3 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 3 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:Chain 4 Exception: lkj_corr_cholesky_lpdf: Random variable[2] is 0, but must be positive! (in '/tmp/Rtmpg5kVA9/model-b56e6f41d818.stan', line 29, column 2 to column 33)Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.Chain 4 Chain 4 finished in 6.3 seconds.
Chain 1 finished in 6.7 seconds.
Chain 3 finished in 6.7 seconds.
Chain 2 finished in 7.6 seconds.
All 4 chains finished successfully.
Mean chain execution time: 6.8 seconds.
Total execution time: 7.6 seconds.my_summary(vgrp_fit$draws())# A tibble: 52 × 10
variable mean sd q2.5 q5 q25 q50 q75 q95
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ 3.58e+2 4.48 3.48e+2 3.50e+2 3.55e+2 3.58e+2 3.61e+2 3.65e+2
2 gamma[1,1] 5.02e-1 0.0551 3.91e-1 4.14e-1 4.66e-1 5.02e-1 5.39e-1 5.92e-1
3 gamma[2,1] 6.18e-1 0.0154 5.87e-1 5.93e-1 6.08e-1 6.18e-1 6.28e-1 6.43e-1
4 gamma[1,2] 4.35e-2 0.0660 -8.68e-2 -6.03e-2 4.12e-4 4.31e-2 8.54e-2 1.54e-1
5 gamma[2,2] 8.01e-2 0.0190 4.29e-2 4.89e-2 6.79e-2 8.04e-2 9.24e-2 1.11e-1
6 sigma 9.26e-2 0.00476 8.37e-2 8.51e-2 8.94e-2 9.24e-2 9.56e-2 1.01e-1
7 tau[1] 6.70e-2 0.0486 2.96e-3 5.76e-3 2.92e-2 6.04e-2 9.40e-2 1.55e-1
8 tau[2] 2.39e-2 0.0131 1.81e-3 3.96e-3 1.55e-2 2.28e-2 3.07e-2 4.66e-2
9 z[1,1] 4.33e-1 0.809 -1.27e+0 -9.48e-1 -5.01e-2 4.40e-1 9.42e-1 1.75e+0
10 z[2,1] 1.46e-1 0.756 -1.45e+0 -1.18e+0 -3.06e-1 1.65e-1 6.41e-1 1.32e+0
# ℹ 42 more rows
# ℹ 1 more variable: q97.5 <dbl>vgrp_fit$diagnostic_summary()$num_divergent
[1] 0 0 0 0
$num_max_treedepth
[1] 0 0 0 0
$ebfmi
[1] 0.7879778 0.8287131 0.8046688 0.791515310.2 brms
priors <- c(
prior(normal(0, 2.5), class = Intercept),
prior(normal(0, 2.5), class = b),
prior(lkj(2), class = cor), # Ω ~ LKJ(2)
prior(cauchy(0, 2.5), class = sigma)
)
vgrp_fit_brm <- brm(
log(h) ~ log(dbh) * wd + (1 + log(dbh) | sp),
family = gaussian(),
data = allo_df,
prior = priors,
refresh = 0, # don't print update
seed = 1234,
chains = 4, cores = 4,
control = list(
adapt_delta = 0.95, # default is 0.8
max_treedepth = 15 # default is 10
),
backend = "cmdstanr"
)Start samplingRunning MCMC with 4 parallel chains...
Chain 1 finished in 6.2 seconds.
Chain 4 finished in 6.7 seconds.
Chain 3 finished in 6.8 seconds.
Chain 2 finished in 7.3 seconds.
All 4 chains finished successfully.
Mean chain execution time: 6.7 seconds.
Total execution time: 7.3 seconds.Warning: 1 of 4000 (0.0%) transitions ended with a divergence.
See https://mc-stan.org/misc/warnings for details.Divergent transitions! We need to increase adapt_delta and max_treedepth to avoid divergent transitions. We may also need to change the priors. Reparameterization is less straightforward in brms than writing your own Stan code.
summary(vgrp_fit_brm)Warning: There were 1 divergent transitions after warmup. Increasing
adapt_delta above 0.95 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup Family: gaussian
Links: mu = identity
Formula: log(h) ~ log(dbh) * wd + (1 + log(dbh) | sp)
Data: allo_df (Number of observations: 200)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Multilevel Hyperparameters:
~sp (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sd(Intercept) 0.09 0.07 0.00 0.25 1.00 1710
sd(logdbh) 0.08 0.03 0.04 0.15 1.00 1533
cor(Intercept,logdbh) 0.15 0.43 -0.72 0.85 1.00 1018
Tail_ESS
sd(Intercept) 1911
sd(logdbh) 2119
cor(Intercept,logdbh) 1332
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.50 0.06 0.39 0.62 1.00 3881 2832
logdbh 0.60 0.03 0.54 0.66 1.00 1778 2069
wd 0.02 0.07 -0.12 0.17 1.00 3192 2064
logdbh:wd -0.05 0.04 -0.13 0.02 1.00 2444 2240
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.09 0.00 0.08 0.10 1.00 5598 3118
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).10.3 lme4
grp_fit_lme <- lme4::lmer(
log(h) ~ log(dbh) * wd + (1 + log(dbh) | sp),
data = allo_df)boundary (singular) fit: see help('isSingular')Not convergent! It’s hard to fit this kind of complex multilevel model with MLE framework like lme4.
11 References
Betancourt, M. (2017). Diagnosing Biased Inference with Divergences.
Gelman, A. et al. Bayesian Data Analysis, Third Edition. (Chapman & Hall/CRC, 2013)
devtools::session_info()─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.5.1 (2025-06-13)
os Ubuntu 24.04.3 LTS
system aarch64, linux-gnu
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz Etc/UTC
date 2025-11-03
pandoc 3.8.1 @ /usr/bin/ (via rmarkdown)
quarto 1.7.32 @ /usr/local/bin/quarto
─ Packages ───────────────────────────────────────────────────────────────────
! package * version date (UTC) lib source
P abind 1.4-8 2024-09-12 [?] RSPM
P backports 1.5.0 2024-05-23 [?] RSPM
P bayesplot * 1.14.0 2025-08-31 [?] RSPM (R 4.5.1)
P boot 1.3-32 2025-08-29 [?] RSPM (R 4.5.1)
P bridgesampling 1.1-2 2021-04-16 [?] RSPM (R 4.5.1)
P brms * 2.23.0 2025-09-09 [?] RSPM (R 4.5.1)
P Brobdingnag 1.2-9 2022-10-19 [?] RSPM (R 4.5.1)
P cachem 1.1.0 2024-05-16 [?] RSPM
P checkmate 2.3.3 2025-08-18 [?] RSPM
P cli 3.6.5 2025-04-23 [?] RSPM
P cmdstanr * 0.9.0 2025-03-30 [?] repository (https://github.com/stan-dev/cmdstanr@da99e2b)
P coda 0.19-4.1 2024-01-31 [?] RSPM (R 4.5.1)
codetools 0.2-20 2024-03-31 [3] CRAN (R 4.5.1)
P data.table 1.17.8 2025-07-10 [?] RSPM
devtools 2.4.6 2025-10-03 [2] RSPM (R 4.5.0)
P digest 0.6.37 2024-08-19 [?] RSPM
P distributional 0.5.0 2024-09-17 [?] RSPM
P dplyr * 1.1.4 2023-11-17 [?] RSPM (R 4.5.1)
ellipsis 0.3.2 2021-04-29 [2] RSPM (R 4.5.0)
P evaluate 1.0.5 2025-08-27 [?] RSPM
P farver 2.1.2 2024-05-13 [?] RSPM (R 4.5.1)
P fastmap 1.2.0 2024-05-15 [?] RSPM
P forcats * 1.0.1 2025-09-25 [?] RSPM (R 4.5.1)
P fs 1.6.6 2025-04-12 [?] RSPM
P generics 0.1.4 2025-05-09 [?] RSPM
P ggplot2 * 4.0.0 2025-09-11 [?] RSPM (R 4.5.1)
P glue 1.8.0 2024-09-30 [?] RSPM
P gridExtra 2.3 2017-09-09 [?] RSPM (R 4.5.1)
P gtable 0.3.6 2024-10-25 [?] RSPM (R 4.5.1)
P hms 1.1.4 2025-10-17 [?] RSPM (R 4.5.1)
P htmltools 0.5.8.1 2024-04-04 [?] RSPM
htmlwidgets 1.6.4 2023-12-06 [2] RSPM (R 4.5.0)
P inline 0.3.21 2025-01-09 [?] RSPM (R 4.5.1)
P isoband 0.2.7 2022-12-20 [?] RSPM (R 4.5.1)
P jsonlite 2.0.0 2025-03-27 [?] RSPM
P kableExtra 1.4.0 2024-01-24 [?] RSPM (R 4.5.1)
P knitr * 1.50 2025-03-16 [?] RSPM
P labeling 0.4.3 2023-08-29 [?] RSPM (R 4.5.1)
lattice 0.22-7 2025-04-02 [3] CRAN (R 4.5.1)
P lifecycle 1.0.4 2023-11-07 [?] RSPM
P lme4 1.1-37 2025-03-26 [?] RSPM (R 4.5.1)
P loo 2.8.0 2024-07-03 [?] RSPM (R 4.5.1)
P lubridate * 1.9.4 2024-12-08 [?] RSPM (R 4.5.1)
P magrittr 2.0.4 2025-09-12 [?] RSPM
MASS 7.3-65 2025-02-28 [3] CRAN (R 4.5.1)
P Matrix 1.7-4 2025-08-28 [?] RSPM (R 4.5.1)
P matrixStats 1.5.0 2025-01-07 [?] RSPM
P memoise 2.0.1 2021-11-26 [?] RSPM
mgcv 1.9-3 2025-04-04 [3] CRAN (R 4.5.1)
P minqa 1.2.8 2024-08-17 [?] RSPM (R 4.5.1)
P mvtnorm 1.3-3 2025-01-10 [?] RSPM (R 4.5.1)
nlme 3.1-168 2025-03-31 [3] CRAN (R 4.5.1)
P nloptr 2.2.1 2025-03-17 [?] RSPM (R 4.5.1)
P patchwork * 1.3.2 2025-08-25 [?] RSPM (R 4.5.1)
P pillar 1.11.1 2025-09-17 [?] RSPM
P pkgbuild 1.4.8 2025-05-26 [?] RSPM
P pkgconfig 2.0.3 2019-09-22 [?] RSPM
pkgload 1.4.1 2025-09-23 [2] RSPM (R 4.5.0)
P plyr 1.8.9 2023-10-02 [?] RSPM (R 4.5.1)
P posterior 1.6.1 2025-02-27 [?] RSPM
P processx 3.8.6 2025-02-21 [?] RSPM
P ps 1.9.1 2025-04-12 [?] RSPM
P purrr * 1.1.0 2025-07-10 [?] RSPM
P QuickJSR 1.8.1 2025-09-20 [?] RSPM (R 4.5.1)
P R6 2.6.1 2025-02-15 [?] RSPM
P rbibutils 2.3 2024-10-04 [?] RSPM (R 4.5.1)
P RColorBrewer 1.1-3 2022-04-03 [?] RSPM (R 4.5.1)
P Rcpp * 1.1.0 2025-07-02 [?] RSPM
P RcppParallel 5.1.11-1 2025-08-27 [?] RSPM (R 4.5.1)
P Rdpack 2.6.4 2025-04-09 [?] RSPM (R 4.5.1)
P readr * 2.1.5 2024-01-10 [?] RSPM (R 4.5.1)
P reformulas 0.4.1 2025-04-30 [?] RSPM (R 4.5.1)
remotes 2.5.0 2024-03-17 [2] RSPM (R 4.5.0)
renv 1.1.5 2025-07-24 [1] RSPM (R 4.5.1)
P reshape2 1.4.4 2020-04-09 [?] RSPM (R 4.5.1)
P rlang 1.1.6 2025-04-11 [?] RSPM
P rmarkdown 2.30 2025-09-28 [?] RSPM
P rstan 2.32.7 2025-03-10 [?] RSPM (R 4.5.1)
P rstantools 2.5.0 2025-09-01 [?] RSPM (R 4.5.1)
P rstudioapi 0.17.1 2024-10-22 [?] RSPM
P S7 0.2.0 2024-11-07 [?] RSPM (R 4.5.1)
P scales 1.4.0 2025-04-24 [?] RSPM (R 4.5.1)
sessioninfo 1.2.3 2025-02-05 [2] RSPM (R 4.5.0)
P StanHeaders 2.32.10 2024-07-15 [?] RSPM (R 4.5.1)
P stringi 1.8.7 2025-03-27 [?] RSPM
P stringr * 1.5.2 2025-09-08 [?] RSPM
P svglite 2.2.1 2025-05-12 [?] RSPM (R 4.5.1)
P systemfonts 1.3.1 2025-10-01 [?] RSPM
P tensorA 0.36.2.1 2023-12-13 [?] RSPM
P textshaping 1.0.4 2025-10-10 [?] RSPM
P tibble * 3.3.0 2025-06-08 [?] RSPM
P tictoc * 1.2.1 2024-03-18 [?] RSPM (R 4.5.1)
P tidyr * 1.3.1 2024-01-24 [?] RSPM (R 4.5.1)
P tidyselect 1.2.1 2024-03-11 [?] RSPM (R 4.5.1)
P tidyverse * 2.0.0 2023-02-22 [?] RSPM (R 4.5.1)
P timechange 0.3.0 2024-01-18 [?] RSPM (R 4.5.1)
P tzdb 0.5.0 2025-03-15 [?] RSPM (R 4.5.1)
usethis 3.2.1 2025-09-06 [2] RSPM (R 4.5.0)
P utf8 1.2.6 2025-06-08 [?] RSPM
P vctrs 0.6.5 2023-12-01 [?] RSPM
P viridisLite 0.4.2 2023-05-02 [?] RSPM (R 4.5.1)
P withr 3.0.2 2024-10-28 [?] RSPM
P xfun 0.54 2025-10-30 [?] RSPM
P xml2 1.4.1 2025-10-27 [?] RSPM (R 4.5.1)
P yaml 2.3.10 2024-07-26 [?] RSPM
[1] /workspaces/afec-bayes-2025/renv/library/linux-ubuntu-noble/R-4.5/aarch64-unknown-linux-gnu
[2] /usr/local/lib/R/site-library
[3] /usr/local/lib/R/library
* ── Packages attached to the search path.
P ── Loaded and on-disk path mismatch.
──────────────────────────────────────────────────────────────────────────────