BayesianRegressionModels — Example Catalog

Binomial regression with logit link (default). Family: binomial, link="logit".

Binomial regression with probit link. Family: binomial, link="probit".

Binomial regression with complementary log-log link. Family: binomial, link="cloglog".

Linear regression on synthetic data; demonstrates blackjax, numpyro and nutpie backends. Family: gaussian.

Intercept-only beta regression on synthetic probability data. Family: beta.

Beta regression of probability on coin dirt differential. Family: beta.

Intercept-only beta regression on baseball batting averages. Family: beta.

Beta regression of batting average on prior-year batting average. Family: beta.

Categorical regression on three Gaussian classes. Family: categorical.

Categorical regression predicting iris species from morphological measurements. Family: categorical.

Categorical regression of alligator food choice on body length and sex. Family: categorical.

Spline basis regression of bloom day-of-year on year, with explicit intercept term in the basis expansion. Family: gaussian.

Spline basis regression of bloom day-of-year on year, with model intercept absorbed (not explicit in the basis). Family: gaussian.

Circular regression of periwinkle direction on distance using von Mises distribution. Family: vonmises.

Gaussian regression of circular direction data; used as a misspecified comparison model against the von Mises. Family: gaussian.

Poisson regression of roach counts with log-exposure offset. Family: poisson.

Negative binomial regression of roach counts with log-exposure offset; accounts for overdispersion relative to Poisson. Family: negativebinomial.

Mixed-effects model of pig weight over time with random slope and intercept per pig. Family: gaussian.

Gamma regression with a single (constant) shape parameter alpha. Family: gamma.

Distributional gamma regression: compound formula models both mu and alpha as functions of x. Family: gamma.

Distributional gamma regression on hourly bike count data; both mu and alpha modeled with spline basis functions. Family: gamma.

Naive logistic regression ignoring group structure entirely. Family: bernoulli.

Fixed effects logistic regression with a dummy indicator per group. Family: bernoulli.

Multilevel logistic regression with random group intercepts. Family: bernoulli.

Mundlak machine: random intercepts plus group-mean of x to separate within/between effects. Family: bernoulli.

Non-hierarchical binomial model: separate intercept per player. Family: binomial.

Hierarchical binomial model with random player intercepts. Family: binomial.

Hilbert-space Gaussian process approximation on 1D synthetic data. Family: gaussian.

HSGP with centered parameterization on 1D synthetic data. Family: gaussian.

HSGP on GAM data with group-specific processes sharing covariance parameters. Family: gaussian.

HSGP on GAM data with group-specific processes and independent covariance parameters. Family: gaussian.

Isotropic 2D Hilbert-space Gaussian process on synthetic data. Family: gaussian.

2D HSGP with group-specific processes sharing covariance parameters. Family: gaussian.

2D HSGP with group-specific processes and independent covariance parameters. Family: gaussian.

Anisotropic 2D HSGP allowing different length scales per dimension. Family: gaussian.

Poisson regression combining fixed effects, random site intercepts, and a 2D spatial HSGP term varying by year. Family: poisson.

Logistic regression of Clinton vote on party ID and party-age interaction. Family: bernoulli.

Multilevel binomial regression with crossed random effects for state, ethnicity, education and interactions; used for post-stratification. Family: binomial.

Linear logistic regression of income on sex, race, age, and weekly hours worked. Family: bernoulli.

Quadratic logistic regression adding squared age and hours terms. Family: bernoulli.

Cubic logistic regression adding cubic age and hours terms. Family: bernoulli.

Negative binomial regression of student absences on program type and math score. Family: negativebinomial.

Negative binomial regression adding program-by-math interaction term. Family: negativebinomial.

Cumulative ordinal model with only threshold parameters, no predictors. Family: cumulative.

Cumulative ordinal model with action, intention, contact and their interactions. Family: cumulative.

Sequential ratio (sratio) ordinal model predicting years at company from total working years. Family: sratio.

Explicit (raw) quadratic polynomial on projectile data. Family: gaussian.

Orthogonal polynomial degree 2 on projectile data. Family: gaussian.

Linear regression of fuel efficiency on horsepower; baseline model. Family: gaussian.

Orthogonal polynomial degree 2 regression of mpg on horsepower. Family: gaussian.

Orthogonal polynomial regression template for degree ∈ 1:9; used for model comparison. Family: gaussian.

Zero-inflated Poisson model of fish catch count. Family: zeroinflatedpoisson.

Logistic regression of Titanic survival on passenger class, sex and age with interactions. Family: bernoulli.

Linear regression of fuel efficiency on horsepower, weight, cylinders and gear. Family: gaussian.

Negative binomial regression of student absences; dataset requires Stata loading (see :negative_binomial). Placeholder DataFrame() used here. Family: negativebinomial.

Logistic regression of "certified fresh" rating on film length and genre style. Family: bernoulli.

Distributional gamma regression on synthetic data; same model as :distributionalvaralpha. Family: gamma.

Logistic regression of well-switching on distance, arsenic level and education. Family: bernoulli.

Logistic regression adding education interaction terms with distance and arsenic. Family: bernoulli.

Polynomial regression using inline transformation I(t**2) to recover gravitational constant from falling ball data. Family: gaussian.

Same as :polyfallingexplicit but using Bambi's alternative curly-brace inline transformation syntax. Family: gaussian.

Polynomial regression using a pre-computed tsquared column rather than inline transformation. Family: gaussian.

Explicit quadratic polynomial regression on projectile data with both linear and quadratic terms. Family: gaussian.

Quadratic polynomial regression using Bambi's poly() function with raw (non-orthogonal) polynomials. Family: gaussian.

Planet-specific quadratic coefficient interacted with Planet to recover different gravitational constants per planet. Family: gaussian.

Mixed-effects model predicting FVC with random slopes per patient; demonstrates out-of-sample prediction for new groups. Family: gaussian.

Spline regression targeting the 10th percentile of BMI. Family: asymmetriclaplace, kappa=0.1.

Spline regression targeting the 50th percentile (median) of BMI. Family: asymmetriclaplace, kappa=0.5.

Spline regression targeting the 90th percentile of BMI. Family: asymmetriclaplace, kappa=0.9.

Spline regression using Gaussian family; used as comparison model to quantile approach. Family: gaussian.

Complete pooling model: single intercept per floor type, ignoring county structure. Family: gaussian.

No-pooling model: separate intercept for each county-floor combination. Family: gaussian.

Partial pooling with random county intercepts, no floor predictor. Family: gaussian.

Partial pooling with floor as fixed effect and random county intercepts. Family: gaussian.

Partial pooling with random slope and intercept for floor by county. Family: gaussian.

Partial pooling with floor and county-level uranium predictor (log_u). Family: gaussian.

Reaction rate modeled by race × object interaction with crossed random effects per subject, without target random effects. Family: gaussian.

Reaction rate modeled by race × object interaction with crossed random effects for both subject and target. Family: gaussian.

Binary shoot/don't-shoot decision modeled by race × object interaction with crossed random effects for subject and target. Family: bernoulli.

Mixed-effects model of reaction time on days of sleep deprivation with random slope and intercept per subject. Family: gaussian.

Simple model with condition as fixed effect and random participant intercepts. Family: gaussian.

Full model with condition, age, gender, random participant intercepts, and crossed random effects for study and stimulus. Family: gaussian.

Weibull survival model on synthetic data with treatment and age effects. Family: weibull.

Weibull survival model of employee tenure with fixed effects only. Family: weibull.

Weibull survival model of employee tenure with random intercepts for field. Family: weibull.

Bernoulli survival model with complementary log-log link on synthetic discrete-time data. Family: bernoulli, link="cloglog".

Binomial discrete-time survival model on Swedish child mortality data. Family: binomial, link="cloglog".

Binomial discrete-time survival with spline baseline hazard. Family: binomial, link="cloglog".

Poisson formulation of discrete-time survival with log-offset for person-time. Family: poisson.

Exponential survival model with only an intercept (constant hazard). Family: exponential, link="log".

Exponential survival model with cat color as a predictor. Family: exponential, link="log".

Standard linear regression on data with outliers; used as baseline. Family: gaussian.

Robust linear regression using Student's t family, which down-weights outliers. Family: t.

Wald (inverse Gaussian) regression of claim costs on age category (as categorical), gender and area. Family: wald, link="log".

Gamma regression of claim costs on age category (numeric), gender and area. Family: gamma, link="log".

ZIP mu component: models the expected count given non-zero. Family: zeroinflatedpoisson.

ZIP psi component: models the probability of excess zeros. Family: zeroinflatedpoisson.

Hurdle Poisson mu component: models the expected count given positive. Family: hurdle_poisson.

Hurdle Poisson psi component: models the probability of a non-zero outcome. Family: hurdle_poisson.

Compact syntax: one formula for both responses; correlated random effects via shared labels p and q; set_rescor(TRUE) adds residual correlations.

Equivalent explicit syntax with separate bf() per response; same model as :btdata_compact but with individually specified sub-formulas.

Extended model: skew-normal family for tarsus, smooth spline for hatchdate in both sub-models.

Standard binomial GLMM for contagious bovine pleuropneumonia incidence across 4 periods.

Custom beta_binomial2 family; vint passes integer auxiliary data (herd size) to the custom family. Accounts for extra-binomial overdispersion.

Distributional model for both mean and log(sigma) as functions of group; demonstrates that the two groups differ in both location and scale.

Negative binomial model for fishing trip catch counts; no zero-inflation component.

Zero-inflated Poisson; zi sub-model predicts structural-zero probability from child.

Both mean and log(sigma) modeled as smooth splines; shared random group effects.

Monotonic main effect of ordered income; mo() constrains the response to be monotonically increasing or decreasing across ordered categories.

Comparison: standard numeric predictor, ignoring the ordinal structure of income.

Comparison: unordered factor coding, ignores both ordinal structure and monotonicity.

Interaction of monotonic income effect with age group.

Adds city-level random slopes for the monotonic income effect.

Fit to multiply-imputed datasets via brm_multiple; standard regression template on complete cases from each imputed dataset.

Joint model: simultaneously imputes chl while modelling bmi; mi() marks imputed variables; set_rescor(FALSE) removes residual correlations between equations.

Variant: chl measured with known standard error se_chl; propagates measurement uncertainty into the joint imputation model.

Nonlinear exponential model; nl=TRUE enables custom nonlinear predictors; b1 + b2 ~ 1 declares both as intercept-only submodels.

Standard linear model on the same exponential-decay data; baseline comparison.

Nonlinear Weibull loss development curve; ult (ultimate loss) varies by accident year; omega = shape, theta = scale.

Phylogenetic random effect only; A is the phylogenetic covariance matrix.

Adds a species-level residual random effect on top of the phylogenetic effect; (1|phylo) captures species-specific deviations not explained by phylogeny.

Phylogenetic effect plus within-species random effect for 5 repeated observations per species.

Meta-analytic model; known SE derived from N; se() passes known measurement error; random study intercept captures between-species heterogeneity.

Poisson GLMM with phylogenetic random effect; models count phenotype data.

Intercept-only multivariate model; establishes baseline residual correlation between tarsus length and back coloration.

Full model; p/q labels allow correlated random effects across responses. Sex and hatch date as fixed effects; foster nest and dam as random effects.

Adds sex × hatchdate interaction to the full model.

Varying intercepts AND slopes; café-level covariance between them. The afternoon effect on wait time differs across cafés.

B is a B-spline basis matrix (each column is a basis function); constructed via bs(year, knots=iknots) in R with ~17 basis columns.

Intercept-only binomial; baseline probability of pulling left across all actors and treatments.

Actor-indexed intercepts plus treatment effects; no pooling across actors.

Multilevel model with partial pooling across actors and blocks. Adaptive regularization via hierarchical priors on actor and block intercepts.

Varying slopes; each actor and block has its own treatment-effect vector.

Binomial likelihood for estimating p(water) from a globe-tossing experiment. family = binomial(link = "identity"); 0 + Intercept fixes the link to identity, estimating p directly.

mid = marriage as index; conditioning on marriage (a collider) induces a spurious negative association between age and happiness.

No conditioning on marriage; age has no direct effect on happiness in this causal model.

Distributional model: both mean and log(sigma) have their own sub-models. sigma grows with x, producing heteroscedastic residuals.

Intercept-only Gaussian regression on adult heights (age >= 18, n=352). Used as a prior predictive / baseline model.

Simple linear regression of height on mean-centered weight in adults (age >= 18). weight_c = weight - mean(weight).

Quadratic polynomial regression of height on standardized weight using the full dataset (all ages). weight_s = standardized weight, weight_s2 = weight_s^2.

Cubic polynomial regression of height on standardized weight using the full dataset (all ages). weight_s3 = weight_s^3.

Nonlinear exponential growth model for height on mean-centered weight in adults (age >= 18). weight_c = weight - mean(weight).

Lognormal nonlinear model treating body weight as a scaled cylinder. w = weight, h = height. Models body weight as a scaled cylinder volume.

Intercept-only Poisson (log link); baseline tool count model.

Contact-index intercepts and slopes on log-population. cid = contact level (high/low).

Scientific power-law model with identity link; a, b vary by contact level cid.

Power-law model where a has a Gaussian process prior over geographic coordinates (latitude and longitude). Requires Kline2 dataset.

AR(1) autoregressive model for hare population using the full dataset (missings present in lag variables).

Equivalent lag-1 regression; Hare_1 = lag of Hare. Uses dropmissing dataset (d2, n=20).

Treating lagged Hare as partially observed/missing. Uses the full dataset with first-row Hare_1 as missing.

Bivariate VAR(1); lognormal family. Lynx feeds back on Hare. Uses dropmissing dataset (d2, n=20).

Positive association between neocortex percentage and kcal per gram before controlling for body mass. Uses complete cases (n=12) with non-missing neocortex.perc.

Negative association; larger species produce less calorie-dense milk. Uses complete cases (n=12) with non-missing neocortex.perc.

Both effects revealed simultaneously — the masked association. Uses complete cases (n=12) with non-missing neocortex.perc.

Clade index model with one intercept per primate clade, no overall intercept. Uses complete cases (n=12) with non-missing neocortex.perc.

Joint missing-data model from Ch 15; b = neocortex.perc with missing values treated as partially observed. Uses the full dataset (all 29 species, including those with missing neocortex data).

Do moralizing gods predict / follow the emergence of writing? Rows with missing moralizing_gods or writing are dropped.

poisson(link = identity); n_panda = success; phi = asymptotic rate, k = growth rate, theta = shape; age_s = standardized age.

Baseline model with no treatment effect. Growth proportional to initial height.

Post-treatment bias: treatment effect absorbed by conditioning on fungus (a mediator). Including fungus blocks the path from treatment to h1.

Causal model: treatment effect recovered by excluding fungus (the mediator). This gives an unbiased estimate of the total effect of treatment on growth.

Brain size ~ body mass + group size; phylogeny ignored. All variables log-standardized.

Same as OLS model plus phylogenetic correlation matrix R passed via data2. R is built from a phylogenetic tree, e.g. via ape::vcv.phylo().

No-pooling model; one parameter per tank. Each tank gets its own intercept with no sharing of information.

Partial pooling; adaptive regularization via multilevel model. Tank-level intercepts share a common hyperprior.

Pooled model ignoring continent. log_gdp_std = log GDP standardized to mean, rugged_std_c = ruggedness standardized and centered.

Continent intercepts via index cid; same ruggedness slope for Africa and Other. cid = "Africa" or "Other".

Continent-specific intercepts AND slopes; ruggedness hurts GDP outside Africa but helps within Africa.

Intercept-only cumulative ordinal model; 6 thresholds estimated for the 7-point response scale.

Action, contact, and intention effects (with interactions) on moral acceptability ratings.

Monotonic effect of ordered education; edu_new = integer-coded edu (1–8). mo() constrains the effect to be monotonically increasing.

Additive model; water_cent and shade_cent are mean-centered. No interaction term.

Interaction model: the effect of water on blooms depends on shade level.

gid = gender index; apparent gender gap in admission rates before conditioning on department.

Gender and department effects; gender gap disappears after conditioning on department (Simpson's paradox).

Custom beta-binomial family from Ch 12; vint passes integer auxiliary data (applications count).

Simple regression of standardized divorce rate on standardized median age at marriage. D = standardized Divorce, A = standardized MedianAgeMarriage.

Simple regression of standardized divorce rate on standardized marriage rate. D = standardized Divorce, M = standardized Marriage.

Multiple regression; M's effect vanishes after conditioning on A, revealing a spurious association between marriage rate and divorce.

Mediator regression: age at marriage predicts marriage rate. Part of the causal model where A -> M -> D and A -> D.

Multivariate brms model; simultaneous causal model for D and M with residual correlation disabled via set_rescor(FALSE).

Measurement error model from Ch 15; D_sd is the known standard error of the divorce estimate (approximated as 0.1 * |D| + 0.01). D_obs is the noisy observation.

Hierarchical binomial model; players nested within positions; estimates batting average with partial pooling.

Separate group means for calcium vs. placebo; Gaussian (or student-t for robustness). Outcome is blood-pressure decrease in mm Hg.

Gaussian model with missing values; brms uses listwise deletion by default. 50 MCAR missing observations out of 500.

Left/right censored Gaussian model; cen1 ∈ {"left", "none", "right"}; censored at 85 and 115.

Interval censored Gaussian model; y3 = upper bound of each interval; values rounded to nearest 5.

ANOVA as hierarchical model; partial pooling across companion types for longevity.

Pooled / no-group baseline model for longevity; comparison model with no companion-type structure.

Robust ANOVA with group-specific means AND sigmas; student-t family for outlier robustness.

ANCOVA controlling for thorax length; thorax_c = centered Thorax; common slope across companion groups.

ANHECOVA with random slope for thorax length per companion group; allows group-specific relationship between body size and longevity.

Multiple regression of SAT score on spending and percent taking SAT; percent taking is a suppressor variable.

Adds complementary proportion (propnottakez = 1 - prcnttake); explicit centering via 0 + Intercept.

Adds interactionz = spendz * prcnttakez (standardized product); student-t for robustness.

Poisson log-linear model with main effects and interaction as random effects; models 4×4 contingency table.

Binomial version of the 4×4 hair-eye contingency table; trials(264) = total N; same random-effects structure as Poisson variant.

Cumulative ordinal probit with standardized assets as predictor; Assets_s = (Assets - mean) / sd.

Hierarchical model with varying intercepts and slopes per subject; || = uncorrelated random effects; z-standardized predictors.

Simple linear regression with z-standardized predictors; student-t family for robustness. Uses n=300 adults.

Logistic regression predicting sex from standardized weight; single metric predictor.

Logistic regression with two metric predictors; collinear but jointly informative for predicting sex.

Robust nonlinear logistic model; a = probability of guessing (mixture weight); identity link; estimated via nl = TRUE.

Measurement-error model via | se(); quadratic trend in family size; state-level random slopes. Accounts for known standard error in the outcome.

Cumulative ordinal regression; Year and Length standardized to Year_s and Length_s.

Intercept-only ordinal probit model; thresholds estimated with cumulative family and no predictors.

Ordinal model with discrimination parameter; disc controls category spacing and varies with X; lf() submodel.

Heteroscedastic Gaussian comparison model; sigma varies with X; alternative to ordinal family.

Separate binomial recall probability per condition; no intercept centering. Six conditions with distinct true probabilities.

Pooled single intercept model for memory recall; baseline / comparison model ignoring condition differences.

Two-way random effects ANOVA with interaction; Gaussian family; partial pooling across positions and organizations.

Robust variant of salary ANOVA with heterogeneous residual variance per cell (Pos:Org); student-t family.

Categorical/softmax regression with separate intercept per category; categorical family.

Intercepts-only softmax model; estimates baseline probability per category with no predictors.

Sequential ordinal regression with category-specific slopes via cs(); sratio family. First conditional logistic dataset.

Sequential ordinal regression with category-specific slopes via cs(); sratio family. Second conditional logistic dataset.

Split-plot ANOVA with tillage × fertilizer interaction plus blocking on Field.

Split-plot ANOVA without explicit field blocking; comparison model to assess the contribution of the field random effect.

Hierarchical Bernoulli model with partial pooling across 28 healers; logit link.

Single-group Gaussian model on IQ scores; prior-predictive / baseline before conditioning on group.

Two groups with group-specific mean AND sigma; student-t family for robustness against outliers.

Intercept-only Bernoulli model; identity link estimates p directly from 50 trials with 15 successes.

Separate Bernoulli intercept per mint; no intercept centering. Mint 1: 6/8 heads, mint 2: 2/7.

Poisson GLMM; interaction of standardised baseline × treatment; random patient intercept.

Same Poisson GLMM with upper truncation at 104 (the maximum possible count per visit).

Simplified main-effects Poisson; omits baseline covariate; comparison model.

Poisson model for fishing trip catch counts; no zero-inflation component.

Zero-inflated Poisson; zi sub-model predicts structural-zero probability from child.

Gaussian with group-specific sigma; the sub-formula sigma ~ 0 + x models log(sigma) as a function of x.

Bayesian quantile regression at the 25th percentile using an asymmetric Laplace distribution.

Sequential ratio ordinal model; cs() = category-specific coefficient for treat; sratio("logit") family estimates P(Y > k | Y ≥ k) at each threshold.

Lognormal survival model; cens() handles right-censored observations; random intercept per patient accounts for within-patient correlation between two kidneys.

Weibull loss development curve; ult (ultimate loss) varies by accident year; omega = shape, theta = scale; nl = TRUE enables the nonlinear formula.

Equal-weight multiple membership; student outcome depends equally on two schools.

Weighted multiple membership; w1 + w2 = 1 for each student; weights the school random effects by the time spent at each school.

Tensor-product spline of area × construction year; random intercept per district.

Distributional model: both mean and log(sigma) are smooth functions; ID1 label allows correlation between the two sets of district random effects.

Population-level hierarchical regression of the PVL-Delta learning rate (α ∈ [0,1]) on clinical group (healthy controls, heroin, amphetamine) in the Iowa Gambling Task. (1|subjID) models within-subject random intercepts across trials. Transformed via logistic link; LogitNormal prior on regression coefficients.

Population-level hierarchical regression of prospect-theoretic reward sensitivity (A ∈ [0,1], power parameter in utility transformation u = |r|^A) on clinical group in the Iowa Gambling Task. LogitNormal prior; logistic link.

Population-level hierarchical regression of loss aversion weight (w ∈ [0,∞], scales negative prospect values as -w·|loss|^A) on clinical group in the Iowa Gambling Task. Exponential link; LogNormal prior on regression coefficients.

Population-level hierarchical regression of the softmax inverse temperature (β ∈ [0,∞]; higher → more deterministic deck choice) on clinical group in the Iowa Gambling Task. Exponential link; LogNormal prior on regression coefficients.

Population-level hierarchical regression of the Rescorla-Wagner learning rate (α) on Peters Delusions Inventory (PDI) total score and session number. (1|ID) captures between-subject variability. α is constrained to [0,1] via logistic transformation; regression prior: Normal(0, 0.5) on PDI coefficient.

Note: learning_rate is a latent parameter estimated by running ActionModels on the trial-level JGET data before this regression is applied.

Population-level hierarchical regression of the Rescorla-Wagner action noise (β, response standard deviation) on PDI score and session. (1|ID) captures subject-level variability. β is unconstrained (exponential link); prior: Normal(0, 0.3) on intercept.

Naive lognormal regression for delay distribution estimation (e.g. symptom onset to case notification). Does not account for interval censoring or right truncation of secondary events, producing estimates biased toward shorter delays. The sigma sub-model (sigma ~ 1) estimates the lognormal standard deviation on the log scale. Baseline model; compare with :delay_marginal.

Marginal lognormal model correcting for double interval censoring (primary and secondary event times observed only within daily windows) and right truncation (only cases with secondary event before obs_time are observable). Uses a custom marginal_lognormal brms family implemented via the primarycensored Stan library.

vreal() passes auxiliary real-valued data to the custom likelihood: observation cutoff, primary censoring window, secondary censoring window, and upper delay bound. The likelihood integrates over primary event time within the censoring window.

Hierarchical nowcast jointly modeling three age strata. enw_expectation(...) includes day-of-week and daily random effects per group (.group notation). enw_reference(~1 + (1|age_group)) allows lognormal delay parameters to vary by age group via partial pooling. enw_report(...) accounts for weekend reporting dips shared across groups.

Extends :enw_age_reference with time-varying delay: rw(week, by=age_group) adds a weekly random walk on delay parameters independently per age group, capturing temporal changes in testing and reporting practices (e.g., shifts in test positivity or lab capacity) that affect different age groups differently.

Default epinowcast model for nowcasting right-truncated hospitalization counts. enw_expectation(~0 + (1|day)) models expected final counts via a geometric random walk (daily random intercepts on the log scale). enw_reference(~1, dist='lognormal') fits a static lognormal delay distribution shared across all reference dates. Negative-binomial observation model with overdispersion parameter φ.

Extends the basic model with a non-parametric report date effect: enw_report(~(1|day_of_week)) adds random intercepts for day of the week of the report date, capturing systematic under-reporting on weekends. All other components identical to :enw_basic.

Replaces the parametric lognormal delay with a fully non-parametric discrete-time hazard model. enw_reference(parametric=~0, non_parametric=~0+delay) estimates a separate baseline hazard for each delay value using a Cox proportional hazards formulation. More flexible than lognormal; useful when the delay distribution is multi-modal or irregular.

Mechanistic renewal-process expectation model that jointly nowcasts and estimates the instantaneous reproduction number Rₜ. r=~1+rw(week) specifies a weekly random walk on log(Rₜ). generation_time=gt_pmf provides a fixed discretised generation time PMF (Gamma-derived). observation=~1+(1|day_of_week) models day-of-week ascertainment variation. enw_reference(~1, dist='lognormal') fits a static lognormal reporting delay.

Complete-pooling estimate: single shared batting-average parameter for all 18 players. Binomial / logit. Compare with :bball_nopooling and :bball_partial.

No-pooling model: independent intercept per player, no shrinkage. Binomial / logit.

Hierarchical shrinkage model: player-level random intercepts partially pool estimates toward the population mean. stan_glmer / Binomial / logit. The rstanarm pooling vignette shows this recovers end-of-season averages better than either complete or no pooling.

Canonical stan_glmer example from the rstanarm mixed effects vignette. Herd-level random intercept for disease incidence across 4 time periods. Binomial / logit.

stan_glm with Gamma family and inverse link, modelling blood clotting time as a function of log plasma concentration, thromboplastin lot, and their interaction. Classic Dobson (2002) §9.2 example.

stan_lm with an R² prior on explained variance. Full interaction model of cloud-seeding conditions. Gaussian family, identity link.

stan_polr proportional-odds model of tobacco consumption group as a function of age and alcohol group. Ordered logistic.

stan_betareg for proportion outcome (0–1): fraction of crude oil converted to gasoline as a function of API gravity, vapor pressure, vaporization temperature, and batch factor. Beta / logit (mean), constant precision.

Extends :gasoline_beta by modelling precision (phi) as a function of gravity and pressure. The | separator follows the betareg/rstanarm distributional formula interface. Beta / logit (mean) + log (precision).

Baseline Gaussian model regressing child cognitive score on whether the mother completed high school. Demonstrates default weakly-informative priors.

Regresses child test score on continuous maternal IQ; single continuous predictor. Gaussian / identity.

Additive model combining both maternal predictors. Gaussian / identity.

Full interaction; the slope of maternal IQ on child score is allowed to differ by high-school completion status. Gaussian / identity.

stan_glm with default weakly-informative priors (autoscaled Normal). Demonstrates how rstanarm's prior autoscaling adapts to predictor scale. Gaussian family, identity link.

stan_glm Poisson regression for post-treatment roach trap counts, controlling for pre-treatment baseline and senior-housing status. Trap-exposure days enter as a log offset. Poisson / log.

Same predictors as :roaches_poisson but with the negative binomial family to accommodate overdispersion and excess zeros. stan_glm.nb / log.

stan_aov two-way factorial model with R² prior on variance explained by the protein source × protein level interaction. Gaussian / identity.

stan_lmer refit treating the factorial structure as nested random intercepts, illustrating partial pooling across factor levels. Gaussian / identity.

stan_glm logistic regression of household well-switching on distance to the nearest safe well (in units of 100 metres). Binomial / logit.

Adds arsenic level of the current well as a second predictor for switching. Binomial / logit.

stan_glm binomial model of survey agreement with gender-role statements as a function of education and respondent sex. Grouped binomial with logit link.

Extends the linear education model with a quadratic term to capture diminishing effects at higher education levels. Binomial / logit.

Factor as predictor; treatment contrast by default (reference level = "A").

Cell-means parameterization; no global intercept; each level gets its own mean.

Monotonic effect; F treated as ordered categorical; mo() constrains the effect to be monotonically increasing or decreasing across levels.

Full factorial Gaussian model: main effects of A and B plus A×B interaction.

A nested in B; A effects estimated separately within each level of B.

Bernoulli(logit) regression; 2×2 factorial design with binary outcome pDV.

Naive OLS ignoring measurement error in both predictor and response; biased estimates.

Measurement error in both response and predictor; me() propagates predictor uncertainty; resp_se() propagates response uncertainty and adds residual sigma.

resp_se() passes known measurement error; sigma = FALSE fixes residual SD to 0; random study intercept captures between-study heterogeneity τ.

Varying intercepts and slopes for subjects, uncorrelated; || suppresses the correlation parameter between intercept and slope.

Varying intercepts and correlated slopes for subjects; | estimates the intercept–slope correlation.

Crossed random effects for subjects and items; both have varying intercepts and slopes.

Distributional model: residual sigma modeled with a subject-level random intercept, allowing between-subject variability in residual spread.

Hierarchical / partial pooling; adaptive shrinkage across subjects.

Complete pooling; single intercept ignoring subject identity.

No pooling; one intercept per subject; cell-means parameterization.

Gaussian regression; higher cognitive load → larger pupil dilation.

Bernoulli(logit) regression; set-size effect on recall accuracy.

Lognormal regression; practice effect: RT decreases linearly on log scale.

Lognormal regression; varying intercepts and correlated slopes per subject; Stroop congruence effect estimated per subject.

Two-parameter mixture model for continuous report working memory data (Bays et al. 2009). The response error is decomposed into:

• A von Mises (on-target memory) component with concentration kappa

• A uniform (guessing) component with probability 1 − inv_logit(thetat)

Both thetat (log-odds of remembering) and kappa (von Mises concentration, log scale) are regressed on set size with by-participant random slopes. The 0 + set_size intercept-free parameterisation gives one coefficient per set-size level (treatment of set_size as factor).

bmm model spec: mixture2p(resp_error = "error")

Three-parameter mixture model extending mixture2p with a non-target swap component (Bays et al. 2009). The response error is a mixture of:

• Von Mises centred on the target (memory; thetat, kappa)

• Von Mises centred on a randomly selected non-target (swap; thetant, same kappa)

• Uniform (guessing)

thetant captures log-odds of swapping to a non-target. Requires non-target feature angles (non_target_1–non_target_5) and set_size.

bmm model spec: mixture3p(resp_error="error", nt_features=paste0("non_target_",1:5), set_size="set_size")

Full Interference Measurement Model (Oberauer & Lin 2017) with all four parameters varying across two conditions:

• c — baseline memory strength (item activation)

• a — associative activation (binding strength)

• s — spatial similarity gradient (spread of interference)

• kappa — von Mises precision (log scale)

Each parameter gets one coefficient per condition via 0 + cond. Non-target feature angles (color_item2–color_item5) and spatial distances (dist_item2–dist_item5) enter the likelihood directly (not as predictors in the formula).

bmm model spec: imm(resp_error="resp_error", nt_features=paste0("color_item",2:5), set_size="set_size", nt_distances=paste0("dist_item",2:5), version="full")

Signal Discrimination Model (Oberauer 2023) fit to a three-condition continuous report experiment. The SDM represents memory as a signal on a circular space; the response reflects discrimination between a memory trace and noise.

c (memory strength, real-valued) and kappa (precision, log scale) each receive one coefficient per condition via 0 + cond. The model uses a custom circular-normal likelihood implemented via brms stanvar().

bmm model spec: sdm(resp_error = "y")

Data augmentation approach for estimating total species richness from incomplete sampling. The dataset includes 50 observed species plus 50 pseudospecies (all-zero detection histories); the occupancy probability of pseudospecies reflects the probability that an unobserved species is actually present.

ff_species is flocker's reserved grouping keyword for species in augmented models. Occupancy random intercept by species captures species-level variation in commonness; detection random slopes by species capture variation in detectability.

flocker call: flock(f_occ=~(1|ff_species), f_det=~uc1+ec1+(1+uc1+ec1|ff_species), flocker_data=fd, augmented=TRUE)

Autologistic dynamic occupancy model: colonisation and persistence share predictor structure; a constant logit-scale offset (auto: ~ 1) is added to persistence probability relative to colonisation. Initial occupancy at equilibrium.

multiseason = "autologistic", multi_init = "equilibrium" in flocker. The autologistic formulation is more parsimonious than explicit colex when persistence and colonisation share the same predictors.

flocker call: flock(f_det=~uc1+ec1, f_col=~uc1, f_auto=~1, flocker_data=fd, multiseason="autologistic", multi_init="equilibrium")

Dynamic (multi-season) colonisation-extinction occupancy model with explicit estimation of initial occupancy. Four sub-models:

• occ — initial occupancy at season 1 (log-odds)

• det — detection given presence (log-odds; with event covariate)

• col — colonisation probability (log-odds): unoccupied sites becoming occupied

• ex — extinction probability (log-odds): occupied sites becoming unoccupied

multi_init = "explicit" in flocker: initial occupancy estimated freely.

flocker call: flock(f_occ=~uc1, f_det=~uc1+ec1, f_col=~uc1, f_ex=~uc1, flocker_data=fd, multiseason="colex", multi_init="explicit")

Dynamic colonisation-extinction occupancy model where initial occupancy is derived from the colonisation/extinction equilibrium rather than estimated as a free parameter. No occ sub-model; initial occupancy = col / (col + ex).

multi_init = "equilibrium" in flocker. Reduces parameters and enforces stationarity of the occupancy process.

flocker call: flock(f_det=~uc1+ec1, f_col=~uc1, f_ex=~uc1, flocker_data=fd, multiseason="colex", multi_init="equilibrium")

Single-season multi-species occupancy model with site-level (uc1) and visit-level (ec1) covariates.

• occ sub-model: log-odds of occupancy at a site; uc1 as fixed slope and random slope by species.

• det sub-model: log-odds of detection given presence; both uc1 and the event-varying ec1 with random slopes by species.

Species-level random effects pool information across the 10 species. The model returns a brmsfit object; all brms post-processing tools apply.

flocker call: flock(f_occ = ~ uc1 + (1 + uc1 | species), f_det = ~ uc1 + ec1 + (1 + uc1 + ec1 | species), flocker_data = fd)

Single-season multi-species occupancy model with only site-level (uc1) covariates in both occupancy and detection sub-models (no event-varying covariate). Detection probability is constant across visits within a site, enabling an efficient "rep-constant" parameterisation in flocker.

flocker call: flock(f_occ = ~ uc1 + (1 + uc1 | species), f_det = ~ uc1 + (1 + uc1 | species), flocker_data = fd, rep_constant = TRUE)

N-mixture model for imperfect detection with non-parametric smooth effects. The N-mixture likelihood jointly models the imperfect-detection observation process and the true abundance process.

• obs formula: detection probability sub-model; smooth effects of two detection covariates

• trend formula: latent abundance sub-model; smooth abundance covariate + hierarchical random effect for the categorical habitat factor

• family = nmix(): N-mixture likelihood (Royle 2004)

Fit with: mvgam(y~s(det_cov,k=4)+s(det_cov2,k=4), trend_formula=~s(abund_cov,k=4)+s(abund_fac,bs="re"), family=nmix(), data=nmix_data)

Vector Autoregressive (VAR) model for multivariate Gaussian plankton time series with correlated process errors. Environmental forcing (temperature × month tensor product) acts at the latent process level.

• obs formula: no observation-level predictors

• trend formula: shared tensor-product smooth te(temp, month, k=c(4,4)) plus species-specific deviations (by = trend); -1 removes intercept

• trend_model = VAR(cor = TRUE): VAR with correlated cross-series innovations

The by = trend idiom in mvgam indexes each latent trend series, allowing species-specific environmental response curves.

Fit with: mvgam(y~-1, trend_formula=~te(temp,month,k=c(4,4))+te(temp,month,k=c(4,4),by=trend)-1, trend_model=VAR(cor=TRUE), family=gaussian(), data=plankton_train)

Baseline Poisson GLM capturing inter-annual variation in rodent counts via hierarchical (random-effects) smooths on year. s(year_fac, bs = "re") - 1 is the mgcv notation for a random intercept per year level with no overall intercept. No explicit latent trend process.

Fit with: mvgam(count ~ s(year_fac, bs="re") - 1, family=poisson(), data=portal_train, trend_model="None")

State-space Poisson model separating a smooth observation-level NDVI effect from a latent AR(1) trend process. The trend captures residual temporal autocorrelation not explained by the environmental covariate.

• obs formula: intercept-only observation model

• trend formula: smooth NDVI effect at the latent process level (k = 6 basis functions)

• trend_model = AR(): AR(1) latent trend

Fit with: mvgam(count ~ 1, trend_formula = ~ s(ndvi_ma12, k=6), trend_model=AR(), family=poisson(), data=portal_train)

Multi-series Poisson model where all five species share a single latent AR(1) trend, with series-specific observation offsets and a cyclic seasonal spline at the trend level.

• obs formula: series - 1 gives a species-specific intercept (fixed offset)

• trend formula: cyclic cubic regression spline for seasonality (bs = "cc")

• trend_model = AR() with trend_map forcing all series to share one trend

The shared-trend structure is specified via a trend_map data.frame passed separately to mvgam() mapping each series to the same latent trend index.

Fit with: mvgam(count~series-1, trend_formula=~s(season,bs="cc",k=8), trend_model=AR(), family=poisson(), data=portal_train, trend_map=trend_map)

Beta-distributed salmon survival proportions with an AR(1) latent trend and a time-varying effect of spring upwelling index (CUI.apr).

• obs formula: intercept-only observation model (Beta family)

• trend formula: dynamic(CUI.apr, k=25, scale=FALSE) makes the regression coefficient on upwelling a smooth function of time (random-walk basis with 25 knots)

• trend_model = AR(): AR(1) latent process in addition to the TVP

The dynamic() function in mvgam fits a time-varying parameter (TVP) model where the coefficient evolves as a random walk or Hilbert-space GP over time.

Fit with: mvgam(survival~1, trend_formula=~dynamic(CUI.apr,k=25,scale=FALSE)-1, trend_model=AR(), family=betar(), data=salmon_train)

Binomial GLMM for CBPP incidence using lme4's matrix-response cbind syntax. Contrast with the brms equivalent: incidence | trials(size) ~ period + (1|herd).

glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
      family = binomial, data = cbpp)

Binomial GLMM for contraception use with a random slope for urban within district — allowing the urban/rural gap to vary across districts. use is recoded from Y/N to integer 0/1.

glmer(use ~ age + I(age^2) + livch + urban + (urban | district),
      family = binomial, data = Contraception)

The simplest possible LMM: grand-mean intercept plus a random intercept for each batch. Introductory example in Bates et al. (JSS 2015, §2).

lmer(Yield ~ 1 + (1 | Batch), data = Dyestuff)

LMM with nested random effects. batch/cask expands to (1|batch) + (1|batch:cask); casks are identified locally within each batch, not globally.

lmer(strength ~ (1 | batch/cask), data = Pastes)

LMM with two fully crossed random intercepts (plates and samples). Neither grouping factor is nested within the other.

lmer(diameter ~ 1 + (1 | plate) + (1 | sample), data = Penicillin)

LMM with correlated by-subject random slope and intercept for Days. (Days | Subject) fits one 2×2 random-effects covariance matrix per subject.

lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)

LMM with uncorrelated by-subject random slope and intercept. (Days || Subject) constrains the off-diagonal covariance to zero (diagonal random-effects covariance matrix).

lmer(Reaction ~ Days + (Days || Subject), data = sleepstudy)

Bernoulli GLMM with two fully crossed random intercepts (subjects × items). r2 is recoded to integer 0/1 from Y/N.

glmer(r2 ~ Anger + Gender + btype + situ + (1|id) + (1|item),
      family = binomial, data = VerbAgg)

Bernoulli GLMM for contraception use in Bangladesh. MixedModels.jl uses abs2() (squaring function) inside @formula instead of R's I(age^2). use recoded from Y/N to integer 0/1.

fit(GeneralizedLinearMixedModel,
    @formula(use ~ 1 + age + abs2(age) + livch + urban + (1|district)),
    contraception, Bernoulli())

The canonical introductory LMM in MixedModels.jl: one variance component for batch-to-batch variability, no fixed covariates.

fit(LinearMixedModel, @formula(Yield ~ 1 + (1 | Batch)), dyestuff)

Note: MixedModels.jl's bundled copy uses lowercase yield/batch; the Rdatasets CSV uses Yield/Batch.

LMM with nested random effects: batch/cask expands to (1|batch) + (1|batch:cask). Cask labels are local to each batch.

fit(LinearMixedModel, @formula(strength ~ 1 + (1 | batch/cask)), pastes)

LMM with two fully crossed random intercepts. Neither plate nor sample is nested within the other; both random intercepts appear additively.

fit(LinearMixedModel,
    @formula(diameter ~ 1 + (1 | plate) + (1 | sample)), penicillin)

LMM with correlated by-subject random slope and intercept for Days. Idiomatic MixedModels.jl style writes explicit 1 + intercepts.

fit(LinearMixedModel, @formula(Reaction ~ 1 + Days + (1 + Days | Subject)), sleepstudy)

Same as the correlated model but with the off-diagonal covariance forced to zero. MixedModels.jl uses zerocorr() instead of lme4's (x || g) syntax.

fit(LinearMixedModel,
    @formula(Reaction ~ 1 + Days + zerocorr(1 + Days | Subject)), sleepstudy)

Bernoulli GLMM with two crossed random intercepts (subjects × items). r2 recoded from Y/N to integer 0/1.

fit(GeneralizedLinearMixedModel,
    @formula(r2 ~ 1 + Anger + Gender + btype + situ + (1|id) + (1|item)),
    verbagg, Bernoulli())

First-order random walk latent field on the time index t; the intercept is absorbed into the random walk (hence 0 +). Gaussian family. Equivalent to a Bayesian smoothing spline.

Second-order random walk; smoother than RW1, penalises second differences. Equivalent to cubic spline smoothing. Gaussian family.

Latent AR(1) process over time. Estimates autocorrelation parameter ρ jointly with the observation variance. Gaussian family.

Seasonal latent field with period 12 months (sum-to-zero constraint within each year). Captures multiplicative seasonality on the log scale. Gaussian family.

Intrinsic CAR (Besag) spatial random effect for 506 Boston census tracts. f(id, model = "besag", graph = W) requires an adjacency matrix W (spdep object) passed separately — it is not contained in the tabular data. Gaussian family. Also fitted with "bym" (BYM = ICAR + IID) for comparison.

inla() with family = "gaussian" — the simplest INLA use case, equivalent to lm() but with Bayesian inference via INLA. Demonstrates default weakly-informative priors on regression coefficients.

1D nonparametric regression using a second-order random walk latent field over the continuous predictor range. constr = FALSE drops the sum-to-zero constraint; -1 removes the fixed intercept (absorbed by the random walk). Gaussian family. Compare with spline-based approaches.

f(plate_id, model = "iid") specifies exchangeable random intercepts per assay plate — the INLA equivalent of (1|plate) in lme4. Gaussian family.

Random intercept per subject; subject_id is the numeric version of Subject. Gaussian family. Equivalent to lmer(Reaction ~ Days + (1|Subject)).

Random slopes model: f(subject_id, Days, ...) specifies subject-specific slopes for Days — the INLA way of writing (0 + Days | Subject). Gaussian family.

12-hospital mortality benchmark. f(hospital, model = "iid") with family = "binomial" (trials = n) shrinks small-hospital rates toward the population mean. Demonstrates PC priors on the precision hyperparameter.

inla.surv() wraps right-censored event times for INLA's survival families. Weibull proportional hazards model (family = "weibullsurv"); the shape parameter is estimated as a hyperparameter.

The canonical MCMCglmm animal model: animal random effect with the inverse numerator relationship matrix Ainv (from BTped via inverseA()) separates additive genetic variance from shared rearing-environment variance. ginverse = list(animal = Ainv) must be passed at fitting time. Gaussian family.

Joint model for tarsus length and back colour. trait - 1 gives one intercept per response; us(trait):animal estimates a 2×2 genetic covariance matrix across traits. Allows estimating genetic correlations between morphological traits. ginverse = list(animal = Ainv) required. Gaussian/Gaussian family.

Simple random intercept for rearing nest; estimates shared nest environment contribution to tarsus length. Gaussian family.

Individual chicks vary in both baseline weight and growth rate. ~us(1 + Time):Chick fits a 2×2 unstructured G-matrix per chick capturing the (co)variance of intercepts and slopes. Gaussian family.

Joint longitudinal model for continuous log-bilirubin (Gaussian) and binary ascites (threshold/probit). family = c("gaussian", "threshold"). The residual covariance ~us(trait):units accounts for within-visit correlation between the two outcomes after conditioning on the shared subject random effect.

Intercept-only fixed effects; ~FSfamily random intercept partitions phenoloxidase variance into between-family and within-family components. Gaussian family. Textbook introductory example from the MCMCglmm reference manual.

Ordered threshold (probit) model for horn morphology: polled < scurred < normal. family = "threshold"; residual variance is fixed at 1 for identifiability. ~id captures repeated-individual variation (some animals appear multiple times).

ar1(times + 0 | group) specifies AR(1) autocorrelation across ordered times within each group as a random-effect covariance structure. Gaussian family. Distinguishes glmmTMB from lme4, which does not support non-diagonal random-effect covariance natively.

family = beta_family() with dispformula = ~x models both the mean proportion and the precision (φ) as functions of x. Demonstrates that ignoring heterogeneous dispersion (using dispformula = ~1) leads to worse fit.

glmmTMB with family = poisson and ziformula = ~1 (constant zero-inflation probability). Nest random intercept handles clustering; offset accounts for brood size. Baseline model for owl sibling negotiation calls.

Same formula as :owls_zip but with family = nbinom2 (variance = μ + μ²/k). Handles both overdispersion and excess zeros. AIC typically preferred over the Poisson version.

family = nbinom1 (variance = φμ, linear mean-variance relationship); brood size as a fixed covariate rather than offset. Contrasts NB1 vs NB2 parameterisation.

Hurdle model: ziformula = ~. mirrors the conditional formula in the zero part, combined with family = truncated_nbinom1 for the positive count part. Treats zero-generation and positive-count processes as completely separate.

Zero-inflation probability varies by species (ziformula = ~spp). The conditional mean count also varies by species with a site random intercept. family = nbinom2.

family = genpois (Generalised Poisson) handles both over- and under-dispersion; zero-inflation depends on both species and mining status.

exp(pos + 0 | group) fits a spatial random field with exponential correlation decay. pos must be a numFactor encoding 2D coordinates (x, y). Gaussian family. All observations share one group so the full spatial covariance is estimated as a single 100×100 block. Note: pos = numFactor(x, y) must be computed before calling glmmTMB.

Canonical example from the GLM.jl documentation (and Dobson 1990). Poisson log-linear model for counts in a two-way layout.

glm(@formula(Counts ~ Outcome + Treatment), dobson, Poisson())

OLS regression of national savings rate on demographic and income predictors (50 countries). Classic textbook dataset (Belsley et al. 1980).

lm(@formula(SR ~ Pop15 + Pop75 + DPI + DDPI), LifeCycleSavings)

Binary logistic regression for diabetes diagnosis in Pima Indian women. type is recoded from Yes/No to integer 0/1.

glm(@formula(type ~ npreg + glu + bp + bmi + ped + age), pima, Binomial())

Negative binomial GLM for overdispersed count data (school absent days). Example from the GLM.jl documentation; R equivalent is MASS::glm.nb().

negbin(@formula(Days ~ Eth + Sex + Age + Lrn), quine, LogLink())

Poisson GLM for school absenteeism counts. Overdispersion relative to Poisson motivates the negative-binomial model (:quine_nb). Same formula, different distributional assumption.

glm(@formula(Days ~ Eth + Sex + Age + Lrn), quine, Poisson())

Linear regression of timber volume on tree girth and height; standard OLS via GLM.jl's lm() (equivalent to glm(..., Normal(), IdentityLink())).

lm(@formula(Volume ~ Height + Girth), trees)