Sampling-geometry release. Multiple-membership models with weakly informed members mixed badly under the previous defaults: standard deviations routinely returned bulk ESS in the single digits and R-hat above 2 even on long runs. The causes were the model parameterization, the default scale priors, and the starting values — not the sampler settings. Refit any model whose convergence mattered; estimates from earlier versions are unaffected in expectation but were often not converged.
mm() and hm() random intercepts and slopes are now sampled in
non-centred form (u = sigma * z, z ~ N(0, 1)) instead of the centred
u ~ N(0, sigma^2). The two imply the same posterior; the non-centred form
removes the funnel between an effect and its own standard deviation, which is
what stalled sigma sampling. Correlated (cor = TRUE) and AR blocks are
unchanged in this release.dscaled.gamma(25, 1) — a half-Cauchy(0, 25) whose
mass sat orders of magnitude above any realistic sd and whose tails fed the
funnel above.mm() block the scale additionally undoes membership attenuation.
A member effect reaches the outcome only through the weighted sum
sum_k w_ik u_k, whose standard deviation is sigma * sqrt(sum_k w_ik^2) —
the square root of the Herfindahl index of the weights, or one over the square
root of the effective number of members. A member-level sd therefore lives on a
scale about sqrt(effN) above its footprint on the outcome, and scaling it to
sd(y) alone would shrink it by that factor.prior(..., class = "sd") specifications are unaffected.w() and fn() parameters now start each chain from its own
dispersed draw, taken on the scale of the parameter's prior and held inside
its support. Previously every chain started from one shared value, which makes
R-hat anti-conservative — it assumes overdispersed starting points.fn("smax", kappa = est()) and fn("gmean", p = est())
are kept away from zero, where those functions are singular. kappa now
starts on both sides of zero rather than always at +1, so a min-like and a
max-like posterior are explored symmetrically.inits are applied exactly as given and are never jittered.n.adapt argument (default 1000, was R2jags's 100). JAGS freezes sampler
tuning when adaptation ends, so models with non-conjugate nodes benefit from a
longer adaptation phase.cores > 1 no longer routes through R2jags::jags.parallel, which silently
ignores n.adapt and cannot give each chain its own starting values. Chains
now run through a small parLapply runner and are recombined identically.run = FALSE additionally returns jags.inits for inspection.monitor is now capability-based and defaults to "summary", which retains
estimates and convergence diagnostics without retaining posterior draws."parameters", "random_effects", "weights",
"fitted", "predictive", or "log_lik" are stored once in a compact
posterior::draws_array. "full" selects every family-supported capability.monitor = "jags" is the explicit escape hatch for retaining the complete
raw R2jags object. Logical monitor values are no longer accepted.summary(), tidy(),
glance(), and bmlCompare() expose compact convergence information.This release redesigns the model-specification grammar around the micro-macro framework (weights, aggregation functions, effects) and aligns the package API with brms. All breaking changes are clean breaks — old spellings error with a migration message; there are no deprecation shims.
mm() block now selects an aggregation function fn() — the f in
the framework's theta^{micro,f}(M_it). fn("sum") (default) is the additive
weighted mean; new emergent types read properties of the whole member set:
fn("var", moment = p) — weighted central moments (distributional shape)fn("hhi"), fn("effn"), fn("entropy") — concentration of the weightsfn("threshold", c =, kappa =) — smooth critical-mass terms; the cutpoint
can be estimated: c = est()fn("smax", kappa =) — smooth max/min; kappa = est() makes the
aggregation function itself the estimand (min <-> mean <-> max)fn("gmean", p = est()) — power/CES mean (substitutes vs complements)fn("cov") — weighted covariance of two member attributesfn(~ E((x - E(x))^2)) — user-written reductions with the
expectation operator E() over the normalized member weights; whitelisted
operations transpile to JAGS. Free symbols become estimable shape parameters.w(~ 1/n), w(~ importance, scale = TRUE), w(~ ilogit(b0 + b1*q)). The old fn(w ~ ...) spelling and
the c = normalization flag were removed (c -> scale, matching brms).w() or fn() formula that is not a
data column is a free parameter; the build messages
Estimating parameters: ... so typos surface immediately.name =; single-feature blocks
can be referenced in main-formula interactions: cross-level
(Ax:education) and block x block (Ax:Vx). Bare references error.RE = re(1 + x) — random intercepts + slopes (lme4-style LHS; grouping
implied by id()). Default cor = FALSE (independent; a deliberate
divergence from brms's correlated default); cor = TRUE opts into the
correlated 2x2 with an LKJ-style prior.FE = fe(1 + x, showFE =) — the no-pooling peer (unit dummies + unit-specific
slopes; first unit is the reference).hm() carries structure only: type = and vars = were removed. Fixed
effects of unit-level covariates go in the main formula
(Y ~ gdp + hm(id(cid), RE = re(1 + gdp))), retiring b.hm.bml(formula, data, family = gaussian(), prior = ...) — data is second.gaussian(), bernoulli(), weibull(),
cox(intervals =) (absorbing cox_intervals); strings still accepted.iter, warmup, thin, chains, cores
(replacing n.iter, n.burnin, n.thin, n.chains, parallel).prior() system: class-targeted priors on the natural (SD) scale,
translated to JAGS automatically — prior(normal(0, 5), class = "b"),
classes Intercept, b, sd, cor, sigma, shape, w, fn;
get_prior() lists every settable prior. Raw JAGS strings remain an escape
hatch. The old priors = list argument was removed.b, labeled by term
((Intercept), majority, A_x, A_x:education); internal arrays b,
b.fn.<k> retire b.mm/b.hm. Weight parameters are b.w.<k> (class
"w"); fn shape parameters are fn.<name>.<k> (class "fn") and are
reported in their own "Function parameters" section, never in the coefficient
table. RE variances print as sd(Intercept), sd(x), cor(Intercept,x).as_draws() / as_draws_df() / as_draws_matrix() / as_draws_array()fixef(), ranef(), coef()fitted(), predict(), posterior_predict() (newdata not yet supported)log_lik(), loo(), waic() (gaussian/bernoulli; pointwise log_lik is
monitored in the generated JAGS model)pp_check() (via bayesplot, Suggests), conditional_effects()make_jagscode() / make_jagsdata() (analogues of make_stancode())mm(ar = TRUE) / hm(ar = TRUE)
were removed (migration error): write RE = re(1, ar = TRUE).RE = re(1, ar = year) orders each member's
participations by a numeric time variable and scales the step variance by
the normalized time gaps, u_t ~ N(u_{t-1}, sigma^2 * gap_t / mean(gap)) —
long gaps allow more drift than short ones; equal spacing reproduces the
ar = TRUE model. Requires an intercept-only re(1); time values must be
unique within each member/unit.bml_delta(), bml_realloc(), bml_turnover() — aggregation over
membership spells: the matched (shift-share style) decomposition of a
composition change into reallocation + attribute change + entry/exit, with
the adding-up identity holding exactly. Observed-panel case; the in-model
spell path (estimated weights, reallocation toward high-RE members) is
future work.fn("var"), stacked feature blocks), estimating the
aggregation function itself (fn("smax", kappa = est())), cross-level
interactions via named blocks (Afin:majority), and heterogeneous member
effects (vars(x + x:z) + re(1 + x) in one block). The friendship example
now adjudicates linear-in-means vs. best-friend with loo_compare() and
shows pp_check(); the Boston example uses prior(..., class = "w").fn(w ~ ...) weights spelling; c = normalization flaghm(type =), hm(vars =), hm(name =) (now labels =), hm(showFE =)mm(ar =) / hm(ar =) — moved into re(1, ar = ...)bml(n.iter =, n.burnin =, n.thin =, n.chains =, parallel =, priors =, cox_intervals =); data as last argumentfamily = "Gaussian"-style capitalized strings still resolve, but the
canonical spellings are lowercase / family functionsPackage renamed from rmm to bml (Bayesian Multiple-Membership Multilevel Models)
New syntax for weight functions: The ar parameter has been moved from the fn() specification to the mm() block level for clearer API
fn(w ~ 1/n, c = TRUE, ar = FALSE)fn(w ~ 1/n, c = TRUE) with ar = FALSE at the mm() levelSupport for multiple mmid groups: The package now supports models with multiple membership identifiers, allowing more complex membership structures
Enhanced documentation: Comprehensive documentation added for the coalgov dataset including:
Flexible weight function parameterization: Enhanced support for parameterizing weight functions with covariates and group-specific structures
Per-group random effects: Random effects can now be specified separately for different mmid groups
Improved JAGS code generation: Optimized model string generation for better performance with complex multiple-membership structures
ar parameter moved: Existing code using fn(w ~ ..., ar = TRUE) must be updated to place ar in the mm() block instead
Dataset changes:
schoolnets dataset (including nodedat and edgedat objects)coalgov dataset with enhanced documentation and additional variablesFixed issues with weight function constraints when using multiple mm() blocks
Improved handling of group-level indices in JAGS variable creation
NEWS.md file to track changes to the package