Getting Started with bml

Introduction

The bml package implements Bayesian Multiple-Membership Multilevel Models with Parameterizable Weight Functions. These models are designed for situations where higher-level outcomes depend on the combined influence of multiple lower-level units. Instead of aggregating lower-level data to higher levels (risking aggregation bias) or disaggregating outcomes to lower levels (artificially inflating sample size), bml explicitly models how lower-level units jointly shape higher-level outcomes.

Installation

Install the stable version from CRAN:

install.packages("bml")

You’ll also need JAGS installed on your system. See the installation vignette for details.

Basic Example

Let’s start with a simple example using coalition government data. Each coalition (higher-level unit) is composed of multiple political parties (lower-level units), and we want to model how party characteristics influence coalition outcomes.

library(bml)
data(coalgov)

# Examine the data structure
head(coalgov[, c("gid", "pid", "pname", "n", "finance", "dur_wkb", "event_wkb")])

Understanding the Data Structure

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

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),
      fn = fn(w ~ 1/n, c = TRUE),
      RE = TRUE
    ),
  family = "Weibull",
  data = coalgov,
  seed = 1
)

summary(mod1)

Breaking down the formula:

Model 2: Parameterizing the Weight Function

A key feature of bml is the ability to model aggregation weights as a function of covariates. 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),
      fn = fn(w ~ 1 / (1 + (n - 1) * exp(-(b0 + b1 * pseat))), c = TRUE),
      RE = TRUE
    ),
  family = "Weibull",
  priors = c("b.w ~ dnorm(0,1)"), # weakly informative prior on the weight parameters
  data = coalgov,
  seed = 1,
  monitor = TRUE
)

summary(mod2)

New elements:

Interpretation:

Visualizing Results

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

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

Next Steps

Key Concepts