Getting Started with bml

Understanding the Data Structure

The coalgov dataset is in long format where each row represents a party’s participation in a coalition:

• gid: Government (coalition) identifier
• pid: Party identifier
• n: Number of parties in the coalition
• finance: Party’s financial dependence on member contributions (standardized)
• dur_wkb: Coalition duration in days
• event_wkb: Early termination indicator (1 = terminated early, 0 = censored)

Model 1: Basic Multiple-Membership Model

Our first model examines coalition duration as a function of party-level characteristics, aggregated using equal weights:

mod1 <- bml(
  Surv(dur_wkb, event_wkb) ~ 1 + majority +
    mm(
      id   = id(pid, gid),
      vars = vars(finance),
      w    = w(~ 1/n, scale = TRUE),
      fn   = fn("sum"),
      RE   = TRUE
    ),
  data   = coalgov,
  family = weibull(),
  seed   = 1
)

summary(mod1)

Breaking down the formula:

• Surv(dur_wkb, event_wkb): Survival outcome (duration and event indicator)
• majority: Government-level covariate (binary indicator)
• mm(): Multiple-membership block containing:
  • id = id(pid, gid): Party ID (member) and government ID (group)
  • vars = vars(finance): Party-level attribute
  • w = w(~ 1/n, scale = TRUE): The weights — who counts (equal weights, normalized to sum to 1)
  • fn = fn("sum"): The aggregation function — how contributions combine (the additive weighted mean; this is the default)
  • RE = TRUE: Party-specific random intercepts (shorthand for RE = re(1))

The feature’s coefficient appears in the output as A_finance — an ordinary class-"b" coefficient of the main model.

Model 2: Parameterizing the Weight Function

A key feature of bml is the ability to estimate the aggregation weights from member characteristics. Instead of assuming all parties have equal influence, we can test whether seat share affects a party’s weight in the coalition:

mod2 <- bml(
  Surv(dur_wkb, event_wkb) ~ 1 + majority +
    mm(
      id   = id(pid, gid),
      vars = vars(finance),
      w    = w(~ 1 / (1 + (n - 1) * exp(-(b0 + b1 * pseat))), scale = TRUE),
      fn   = fn("sum"),
      RE   = TRUE
    ),
  data   = coalgov,
  family = weibull(),
  prior  = prior(normal(0, 1), class = "w"),  # weakly informative prior on the weight parameters
  seed   = 1
)

summary(mod2)

New elements:

• The weight formula includes pseat (party seat share). By the one parameter rule, the symbols b0 and b1 are not data columns and therefore become free parameters — the build messages Estimating parameters: b0, b1 (w.1).
• prior = prior(normal(0, 1), class = "w"): A structured prior on the weight parameters, written on the standard-deviation scale and translated to JAGS automatically. Run get_prior(formula, data, family) to list everything that is settable.

Interpretation:

• If \(b_1 \approx 0\): Seat share doesn’t affect weights (equal influence)
• If \(b_1 > 0\): Larger parties have more weight
• If \(b_1 < 0\): Smaller parties have more weight

Model 3: Beyond the sum — emergent features

The aggregation function fn() is not limited to the weighted mean. Emergent features are properties of the whole member set:

mod3 <- bml(
  Surv(dur_wkb, event_wkb) ~ 1 + majority +
    mm(id = id(pid, gid), vars = vars(finance), w = w(~ 1/n), fn = fn("sum"), RE = TRUE) +      # mean:     A_finance
    mm(id = id(pid, gid), vars = vars(finance), w = w(~ 1/n), fn = fn("var")) + # spread:   V_finance
    mm(id = id(pid, gid), w = w(~ pseat, scale = TRUE), fn = fn("hhi"),         # dominance: C_w
       RE = FALSE),
  data   = coalgov,
  family = weibull(),
  seed   = 1
)

Other built-ins include fn("threshold", c = est()) (critical mass with an estimated cutpoint), fn("smax", kappa = est()) (is the aggregation mean-like or max/min-like?), fn("gmean", p = est()) (CES/power mean), and fn("cov"). You can also write your own reduction with the expression DSL: fn(~ E((x - E(x))^2)).

Visualizing Results

# Coefficients (all class-"b" terms, labeled by feature name)
fixef(mod2)

# Diagnostic plot for a weight parameter
monetPlot(mod2, parameter = "b.w.1[1]", label = "Seat share effect")

# MCMC diagnostics
mcmcDiag(mod2, parameters = "b.w.1[1]")

# The posterior toolchain
draws <- as_draws_df(mod2)