Seeds: Random Effect Logistic Regression

This is the classic Seeds example from Volume 1 of the BUGS examples (see also the OpenBUGS write-up). The data come from a seed germination experiment laid out as a 2×2 factorial design across 21 plates. For each plate we record the number of seeds that germinated, r[i], out of the total number of seeds sown, n[i]. Two experimental factors are crossed: the seed type (x1) and the type of root extract used (x2).

The scientific question is how seed type, root extract, and their interaction affect the probability of germination, while acknowledging that plates differ from one another for reasons the covariates do not capture. To handle this extra plate-to-plate variability (over-dispersion relative to a plain binomial model), the model is a random-effects logistic regression: each plate gets its own random intercept b[i], drawn from a common normal distribution whose precision tau is estimated from the data.

This example demonstrates:

• logistic regression for a factorial experiment,
• random effects for extra-binomial variation,
• translating BUGS link-function syntax into Julia-native syntax, and
• starting a model from supplied initial values.

Model

Let $p_i$ be the germination probability on plate $i$. The model is

\[\begin{aligned} b_i &\sim \text{Normal}(0, \tau) \\ \text{logit}(p_i) &= \alpha_0 + \alpha_1 x_{1i} + \alpha_2 x_{2i} + \alpha_{12} x_{1i} x_{2i} + b_i \\ r_i &\sim \text{Binomial}(n_i, p_i) \end{aligned}\]

Here $\alpha_0$ is the baseline log-odds of germination, $\alpha_1$ and $\alpha_2$ are the main effects of seed type and root extract, and $\alpha_{12}$ is their interaction. The plate effects $b_i$ capture variation left unexplained by those factors. The regression coefficients receive vague priors, as does the random-effect precision tau; the derived quantity sigma = 1 / sqrt(tau) is the standard deviation of the plate effects.

Here is the model written with the @bugs macro. Because Julia treats f(x) = ... as a function definition, the BUGS link-function form logit(p[i]) <- ... is written by applying the inverse link (logistic) on the right-hand side. JuliaBUGS can also run the original BUGS program directly; see Migrating from WinBUGS, OpenBUGS, and JAGS for that workflow.

using JuliaBUGS

seeds = @bugs begin
    for i in 1:N
        r[i] ~ dbin(p[i], n[i])
        b[i] ~ dnorm(0.0, tau)
        p[i] = logistic(alpha0 + alpha1 * x1[i] + alpha2 * x2[i] + alpha12 * x1[i] * x2[i] +
                        b[i])
    end
    alpha0 ~ dnorm(0.0, 1.0e-6)
    alpha1 ~ dnorm(0.0, 1.0e-6)
    alpha2 ~ dnorm(0.0, 1.0e-6)
    alpha12 ~ dnorm(0.0, 1.0e-6)
    tau ~ dgamma(0.001, 0.001)
    sigma = 1 / sqrt(tau)
end

BUGSModelDef:
begin
    for i = 1:N
        r[i] ~ dbin(p[i], n[i])
        b[i] ~ dnorm(0.0, tau)
        p[i] = logistic(alpha0 + alpha1 * x1[i] + alpha2 * x2[i] + alpha12 * x1[i] * x2[i] + b[i])
    end
    alpha0 ~ dnorm(0.0, 1.0e-6)
    alpha1 ~ dnorm(0.0, 1.0e-6)
    alpha2 ~ dnorm(0.0, 1.0e-6)
    alpha12 ~ dnorm(0.0, 1.0e-6)
    tau ~ dgamma(0.001, 0.001)
    sigma = 1 / sqrt(tau)
end

Data

The data are supplied as a NamedTuple. r and n are the germinated and total seed counts on each of the N = 21 plates, while x1 and x2 are the (0/1) indicators for seed type and root extract.

data = (
    r = [10, 23, 23, 26, 17, 5, 53, 55, 32, 46, 10, 8, 10, 8, 23, 0, 3, 22, 15, 32, 3],
    n = [39, 62, 81, 51, 39, 6, 74, 72, 51, 79, 13, 16, 30, 28, 45, 4, 12, 41, 30, 51, 7],
    x1 = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
    x2 = [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1],
    N = 21
)

model = seeds(data)

BUGSModel (parameters are in transformed (unconstrained) space, with dimension 26):

  Model parameters:
    alpha2
    b[21], b[20], b[19], b[18], b[17], b[16], b[15], b[14], b[13], b[12], b[11], b[10], b[9], b[8], b[7], b[6], b[5], b[4], b[3], b[2], b[1]
    tau
    alpha12
    alpha1
    alpha0

  Variable sizes and types:
    b: size = (21,), type = Vector{Float64}
    p: size = (21,), type = Vector{Float64}
    n: size = (21,), type = Vector{Int64}
    alpha2: type = Float64
    sigma: type = Float64
    alpha0: type = Float64
    alpha12: type = Float64
    N: type = Int64
    tau: type = Float64
    alpha1: type = Float64
    r: size = (21,), type = Vector{Int64}
    x1: size = (21,), type = Vector{Int64}
    x2: size = (21,), type = Vector{Int64}

All of the classic examples ship with the package under JuliaBUGS.BUGSExamples, bundling the model definition, data, initial values, and published reference results (this one is JuliaBUGS.BUGSExamples.VOLUME_1.seeds).

Initial values

The initial values published with the classic example set the regression coefficients to zero and the random-effect precision to 10:

inits = (alpha0=0.0, alpha1=0.0, alpha2=0.0, alpha12=0.0, tau=10.0)
model = seeds(data, inits)

JuliaBUGS draws the omitted plate effects b from their prior. See Initial Values for partial initialization, array-valued parameters, and the flat vectors accepted by samplers.

Sampling

To draw posterior samples, construct the model with gradient support and run the No-U-Turn sampler:

using AbstractMCMC, AdvancedHMC, ADTypes, Mooncake, FlexiChains
using LogDensityProblems

model = seeds(data; adtype=AutoMooncake(; config=nothing))

n_samples, n_adapts = 2000, 1000
D = LogDensityProblems.dimension(model)
chain = AbstractMCMC.sample(
    model, NUTS(0.8), n_samples;
    chain_type=VNChain, n_adapts=n_adapts,
    init_params=rand(D), discard_initial=n_adapts,
)
summarystats(chain)

Results

The published reference posterior summaries for this example are:

A correctly converged chain's summarystats output should match these values up to Monte Carlo error.

See also: gallery overview, getting-started tutorial, and migration guide.