1What the method is for

Asking whether settlements whose population grew lost more dialect vocabulary — as a posterior distribution over a coefficient.

The general formulation of Bayesian estimation sits in perspective 3 of Seven Perspectives in Equations. This page does not repeat it. It deals with what has to be decided when a hierarchical model is actually fitted. Four things have to be decided: the likelihood (the response distribution and the link function), the prior, the sampler settings, and how convergence and fit are to be checked.

What this study wanted was an answer of the form “the posterior probability that \(\beta > 0\)” for the coefficient \(\beta\) on standardized population growth rate. The response is a rate of vocabulary loss, a proportion lying strictly inside \((0,1)\). Fitting a Gaussian would misspecify the distribution, so a beta distribution with a logit link was used instead (Ferrari & Cribari-Neto 2004). Settlements are nested within middle-school districts, so the intercept for each district was given a random effect.

There are only 25 settlements in the sample. In a small sample the prior matters. That is why the priors were stated explicitly rather than left at their defaults, so that what each choice implies can be said out loud.

2Origin — the principle from the 18th century, the tools from 2017

Two hundred and fifty years passed before the theorem became usable.

What Stan and brms changed was not the principle but the barrier to entry

Before Stan, using a Bayesian hierarchical model meant writing the sampling code yourself, which in practice restricted it to statisticians. brms made it possible to specify a model in almost the same notation as lm(). The implementation discussed on this page is written as the single line brm(loss ~ ... , family = Beta(link="logit")). Applied researchers have only very recently been able to reach for this method directly.

3Where it comes from mathematically

The axioms of probability, Markov chains, and nineteenth-century analytical mechanics.

Bayesian estimation itself rests on Kolmogorov’s (1933) axiomatization of probability. Writing “\(P(\beta > 0) = 0.99\)” is possible because probability is rigorously defined as a measure assigned to sets. That foundation is shared with frequentist procedures; it does not belong to the Bayesian view alone.

The computation comes from a different lineage. The denominator of the posterior is an integral, and in a hierarchical model it has no analytic solution. So instead of finding the distribution itself, the approach switches to drawing samples from it. The mechanism for drawing those samples is a Markov chain (Markov 1906). Using the property that once the present state is fixed the past is irrelevant to the future, a chain is built whose stationary distribution is the posterior, and then allowed to walk.

What is wanted, and what is actually computed
\[ p(\theta \mid D) \;=\; \frac{p(D \mid \theta)\,p(\theta)}{\displaystyle\int p(D\mid\theta)\,p(\theta)\,d\theta} \qquad\longrightarrow\qquad \theta^{(1)}, \theta^{(2)}, \dots, \theta^{(S)} \;\sim\; p(\theta \mid D) \]

The denominator on the left cannot be solved for a hierarchical model. Given a sequence of draws from the posterior, as on the right, the posterior mean, the interval and \(P(\beta>0)\) all come out by counting over the draws. MCMC is a procedure that replaces an integral with samples.

What Stan uses is Hamiltonian mechanics

Naive MCMC is a random walk that gropes its way forward, and it loses efficiency as dimension rises. The Hamiltonian Monte Carlo that Stan adopts treats the log posterior as a landscape of potential energy and rolls a ball across it according to Hamilton’s (1833) equations of motion. Because it uses gradient information, it can travel a long way in one move even in high dimensions. The default is the No-U-Turn Sampler (Hoffman & Gelman 2014), which sets the step size and path length automatically. Nineteenth-century analytical mechanics sits at the core of twenty-first-century Bayesian computation.

4Definitions and assumptions — what is being decided

The likelihood, the prior and the sampler settings are all chosen by the analyst.

The model used in this study
\[ \mathrm{loss}_i \;\sim\; \mathrm{Beta}\bigl(\mu_i \phi,\; (1-\mu_i)\phi\bigr), \qquad \mathrm{logit}(\mu_i) \;=\; \alpha + u_{d(i)} + \beta\, x_i + \gamma_1 z_{1i} + \gamma_2 z_{2i} \]

\(\mu_i\) … the mean loss rate for settlement \(i\), mapped onto the real line by the logit link
\(\phi\) … the precision parameter, which sets the dispersion of the beta distribution
\(u_{d(i)}\) … the random intercept for the middle-school district \(d\) containing settlement \(i\) (five districts)
\(x_i\) … standardized population growth rate. \(z_1, z_2\) … standardized log sample sizes (middle-school students and the older generation)

Why a weakly informative prior

There are three options. A strong prior settles the conclusion in advance. A flat, near-uninformative prior tends to destabilize the sampler in a small sample, and on the logit scale it spreads probability over coefficient values that could not occur in reality. A weakly informative prior lies between the two and supplies only a thin guide of the form “roughly within this range”.

The priors actually used
\[ \beta,\gamma \sim \mathcal{N}(0,\,1), \qquad \alpha \sim \mathcal{N}(0,\,2), \qquad \phi \sim \mathrm{Gamma}(0.1,\,0.1) \]

Placing \(\mathcal{N}(0,1)\) on a coefficient amounts to saying that on the logit scale the effect is probably somewhere between \(-3\) and \(+3\). For standardized predictors that is wide enough, and narrow enough to prevent divergence. Because it is centred on zero, this prior does not work in the direction of making effects look larger.

A prior is not something that has no influence; it is something whose influence is checked

In this study the data carry more information than the prior, so the effect on the final result is slight. But that is something confirmed after the fact, not something that follows automatically from placing a prior. With a sample of n = 25, a narrower prior would certainly pull the estimate towards itself. Priors are stated in the report so that a reader can check this point.

5How it is used in this research

Run four chains, confirm convergence, then count probabilities off the posterior.

  1. Data　The 25 neighbourhood associations for which a loss rate could be obtained, across five middle-school districts. The response loss lies strictly inside \((0,1)\), so no endpoint transformation was needed.
  2. Model specification　brm(loss ~ z_pop_growth + z_log_ss_ms + z_log_ss_gp + (1|district), family = Beta(link="logit")). The log sample sizes of both generations enter as covariates to control for confounding by sampling effort.
  3. Sampler settings　Four chains, 4,000 iterations each, of which 2,000 are warmup, giving 8,000 post-warmup draws in total. adapt_delta = 0.97, max_treedepth = 12, seed = 1. The backend is cmdstanr.
  4. Checking convergence　Inspect \(\hat{R}\), the effective sample size, and the number of divergent transitions. The reported model had \(\hat{R} = 1.00\), bulk-ESS > 3,000, and at most one divergence.
  5. Counting the posterior probability　Count what fraction of the 8,000 draws have \(\beta > 0\). That fraction is \(P(\beta > 0)\).
  6. Checking robustness　A model with added covariates, a model with a spatial random effect, a count-level model propagating within-settlement uncertainty, a leave-one-settlement-out sensitivity analysis, a re-analysis with merged units, and a power simulation were all run alongside.
Model\(\beta\)95% credible interval\(P(\beta>0)\)
base (main model)+0.231[+0.06, +0.40]0.993
cov (area, elevation and built-up land added)+0.229[+0.03, +0.42]0.987
sp (ICAR spatial random effect added)+0.228[+0.04, +0.40]0.992

The coefficient on standardized population growth rate. Figures are quoted from the run log of R/12_brms_loss.R. The main model had \(\hat{R} = 1.00\), Bulk-ESS 5,184, Tail-ESS 5,006. The spatial model is covered on the ICAR spatial random effects page.

A credible interval and a confidence interval differ in where the probability is placed

The interval above runs from the 2.5th to the 97.5th percentile of the posterior. It may be read as “there is a 95% probability that the true value lies in this interval”. The probability is placed on the parameter. A frequentist confidence interval is the reverse: the parameter is treated as a fixed value and the probability is placed on the procedure. The claim is that if the same procedure were repeated, 95% of the intervals it produced would contain the true value; it says nothing about the probability attaching to the one interval in hand.

This difference is not a difference in the principle of estimation. The likelihood is used in Bayesian work too — it is the \(p(D\mid\theta)\) in the expression above, and maximum likelihood is simply one principle for choosing the point that makes that likelihood largest. What separates the Bayesian and frequentist views is not whether the likelihood is used, but whether a parameter is granted a distribution.

Partial pooling by lexical stratum

The vocabulary was divided into three strata — Kansai, Mino and Tokyo — and a growth slope was estimated for each by partial pooling. Kansai \(\beta = 0.233\) [0.073, 0.394], Mino \(\beta = 0.240\) [0.083, 0.403], Tokyo \(\beta = 0.219\) [0.022, 0.402]. The posterior probabilities are 0.997, 0.996 and 0.983 respectively. The slopes are nearly identical because strong shrinkage is at work, which reflects the effects being similar across strata.

6Limits — what the method cannot say

A posterior probability is a probability conditional on the model being right.

Figures produced by a chain that has not converged carry no meaning whatsoever

\(\hat{R}\) is a diagnostic measuring whether several chains are exploring the same distribution; on convergence it approaches 1. Vehtari et al. (2021) showed that the conventional \(\hat{R}\) can miss convergence failure when the distribution is heavy-tailed or when variances differ between chains, and proposed a rank-normalized revision. What current brms reports is that revised version.

In the diversity model of the same study this became an actual problem. Some of the non-spatial models reached \(\hat{R} = 1.24\), and the results were judged not to have converged and were not adopted. In the spatial models max_treedepth was reached every time, and the spatial variance stuck close to zero, giving a funnel-shaped posterior. Both are cases where “a number came out” and “the quantity could be estimated” are different things. The diagnostics have to be read before the posterior summary.

“Posterior probability 0.99” does not mean “robust”

The main result of this study is \(P(\beta>0) = 0.99\), but the author does not present it as a robust conclusion. In the count-level model, which propagates the within-settlement sampling uncertainty, the corresponding interaction is \(\beta = -0.17\) with a 95% credible interval of [\(-0.50\), \(0.18\)], straddling zero. In the leave-one-settlement-out sensitivity analysis, \(\beta\) falls to 0.12 only when the single settlement with just one older-generation respondent is dropped. Merging settlements until each has five or more older-generation respondents and re-analysing makes the effect disappear.

The power simulation explains the discrepancy. Under the actual sampling design, the probability of detecting a true interaction effect of \(-0.15\) or \(-0.30\) is only 12% and 39% respectively (reaching 82% only at \(-0.50\)). The settlement-level analysis returns a strong conclusion because it treats each settlement’s loss rate as a precise quantity measured without error. A high posterior probability shows only how plausible something is under the assumptions the model has made.

The same misuse that afflicts testing also occurs in Bayesian work

Using “the interval excludes zero” as a substitute for significance is doing the same thing under a different name. Where this study wants to claim that an effect is absent, it reports, alongside the interval, how much posterior probability falls inside a region of practical equivalence of \(\pm 0.05\) and \(\pm 0.10\). This is a measure for stating how precisely something can be said, rather than asserting that there is no effect.

7Sources and data

The subject is a three-generation dialect vocabulary survey conducted in settlements running from Maibara in Shiga Prefecture to Sekigahara in Gifu Prefecture. The original survey was carried out in 2015 by Habiro Net, a local citizens’ group, and this research is a secondary analysis of its anonymized data. The analysis is written across scripts R/01 to R/25, in a structure where the original data files are never modified and are read only through a working copy. The paper is in preparation and was not published at the time of writing.

References on the method

  • Bürkner, P.-C. (2017) brms: An R package for Bayesian multilevel models using Stan. Journal of Statistical Software, 80(1), 1–28. doi:10.18637/jss.v080.i01
  • Carpenter, B., Gelman, A., Hoffman, M. D. et al. (2017) Stan: A probabilistic programming language. Journal of Statistical Software, 76(1), 1–32. doi:10.18637/jss.v076.i01
  • Hoffman, M. D. & Gelman, A. (2014) The No-U-Turn Sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15, 1593–1623.
  • Vehtari, A., Gelman, A., Simpson, D., Carpenter, B. & Bürkner, P.-C. (2021) Rank-normalization, folding, and localization: an improved \(\hat{R}\) for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667–718. doi:10.1214/20-BA1221
  • Ferrari, S. & Cribari-Neto, F. (2004) Beta regression for modelling rates and proportions. Journal of Applied Statistics, 31(7), 799–815.

Which perspective this method belongs to

Bayesian estimation grants a distribution to the parameter and updates prior knowledge with the data into a posterior. What it returns is not a point estimate but a distribution, giving an answer directly in the form of “the probability that the effect is positive”. It maps straight onto perspective 3 (Bayesian) in Seven Perspectives in Equations. ICAR spatial random effects, used in the same study, belongs to the same perspective, in that spatial structure is written as a prior.

This page sets out, for readers in the same field, a method Onohara uses in actual research. Every figure is quoted from the run logs and output files of the analysis scripts; no estimated or approximate values are given.
The whole picture of where each method sits is in Seven Perspectives in Equations, and which teaching session falls under which perspective is in Methods & Mathematics.