Epilepsy: Repeated Measures on Poisson Counts

Breslow and Clayton (1993) analysed data originally reported by Thall and Vail (1990) on seizure counts from a randomised trial of anti-convulsant therapy in epilepsy. Fifty-nine patients were followed over four successive clinic visits, and the number of seizures in the interval before each visit was recorded, along with each patient's treatment assignment, baseline seizure count, and age. The question is whether the treatment reduces the seizure rate after adjusting for these covariates.

The model is a Poisson generalised linear mixed model (model III of Breslow and Clayton). The log seizure rate is a linear function of centred covariates — log baseline count, treatment, a treatment-by-baseline interaction, log age, and an indicator for the fourth visit — plus a subject-level random effect and a subject-by-visit random effect that allows for extra-Poisson variation. This example is from Volume 1 of the classic BUGS examples; see also the OpenBUGS documentation page.

Model

For patient $j = 1, \dots, 59$ and visit $k = 1, \dots, 4$,

\[\begin{aligned} y_{jk} &\sim \text{Poisson}(\mu_{jk}) \\ \log \mu_{jk} &= a_0 + \alpha_{\text{Base}} \left( \log(\text{Base}_j / 4) - \overline{\log(\text{Base}/4)} \right) + \alpha_{\text{Trt}} \left( \text{Trt}_j - \overline{\text{Trt}} \right) + \alpha_{\text{BT}} \left( \text{BT}_j - \overline{\text{BT}} \right) \\ &\quad + \alpha_{\text{Age}} \left( \log(\text{Age}_j) - \overline{\log(\text{Age})} \right) + \alpha_{\text{V4}} \left( \text{V4}_k - \overline{\text{V4}} \right) + b1_j + b_{jk} \\ b1_j &\sim \text{Normal}(0, \sigma_{b1}^2), \qquad b_{jk} \sim \text{Normal}(0, \sigma_b^2) \end{aligned}\]

with vague normal priors on the coefficients and vague gamma priors on the random-effect precisions. The original BUGS program uses R-style dotted names such as alpha.Base and tau.b; in Julia these are written with the var"..." syntax, which lets a variable name contain a dot.

using JuliaBUGS

epil = @bugs begin
    for j in 1:N
        for k in 1:T
            mu[j, k] = exp(
                a0 + var"alpha.Base" * (var"log.Base4"[j] - var"log.Base4.bar") +
                var"alpha.Trt" * (Trt[j] - var"Trt.bar") +
                var"alpha.BT" * (BT[j] - var"BT.bar") +
                var"alpha.Age" * (var"log.Age"[j] - var"log.Age.bar") +
                var"alpha.V4" * (V4[k] - var"V4.bar") + b1[j] + b[j, k]
            )
            y[j, k] ~ dpois(mu[j, k])
            b[j, k] ~ dnorm(0.0, var"tau.b")       # subject*visit random effects
        end
        b1[j] ~ dnorm(0.0, var"tau.b1")        # subject random effects
        BT[j] = Trt[j] * var"log.Base4"[j]    # interaction
        var"log.Base4"[j] = log(Base[j] / 4)
        var"log.Age"[j] = log(Age[j])
    end

    # covariate means:
    var"log.Age.bar" = mean(var"log.Age"[:])
    var"Trt.bar" = mean(Trt[:])
    var"BT.bar" = mean(BT[:])
    var"log.Base4.bar" = mean(var"log.Base4"[:])
    var"V4.bar" = mean(V4[:])

    # priors:
    a0 ~ dnorm(0.0, 1.0E-4)
    var"alpha.Base" ~ dnorm(0.0, 1.0E-4)
    var"alpha.Trt" ~ dnorm(0.0, 1.0E-4)
    var"alpha.BT" ~ dnorm(0.0, 1.0E-4)
    var"alpha.Age" ~ dnorm(0.0, 1.0E-4)
    var"alpha.V4" ~ dnorm(0.0, 1.0E-4)
    var"tau.b1" ~ dgamma(1.0E-3, 1.0E-3)
    var"sigma.b1" = 1.0 / sqrt(var"tau.b1")
    var"tau.b" ~ dgamma(1.0E-3, 1.0E-3)
    var"sigma.b" = 1.0 / sqrt(var"tau.b")

    # re-calculate intercept on original scale:
    alpha0 = a0 - var"alpha.Base" * var"log.Base4.bar" - var"alpha.Trt" * var"Trt.bar" -
             var"alpha.BT" * var"BT.bar" - var"alpha.Age" * var"log.Age.bar" -
             var"alpha.V4" * var"V4.bar"
end

BUGSModelDef:
begin
    for j = 1:N
        for k = 1:T
            mu[j, k] = exp(a0 + var"alpha.Base" * (var"log.Base4"[j] - var"log.Base4.bar") + var"alpha.Trt" * (Trt[j] - var"Trt.bar") + var"alpha.BT" * (BT[j] - var"BT.bar") + var"alpha.Age" * (var"log.Age"[j] - var"log.Age.bar") + var"alpha.V4" * (V4[k] - var"V4.bar") + b1[j] + b[j, k])
            y[j, k] ~ dpois(mu[j, k])
            b[j, k] ~ dnorm(0.0, var"tau.b")
        end
        b1[j] ~ dnorm(0.0, var"tau.b1")
        BT[j] = Trt[j] * var"log.Base4"[j]
        var"log.Base4"[j] = log(Base[j] / 4)
        var"log.Age"[j] = log(Age[j])
    end
    var"log.Age.bar" = mean(var"log.Age"[:])
    var"Trt.bar" = mean(Trt[:])
    var"BT.bar" = mean(BT[:])
    var"log.Base4.bar" = mean(var"log.Base4"[:])
    var"V4.bar" = mean(V4[:])
    a0 ~ dnorm(0.0, 0.0001)
    var"alpha.Base" ~ dnorm(0.0, 0.0001)
    var"alpha.Trt" ~ dnorm(0.0, 0.0001)
    var"alpha.BT" ~ dnorm(0.0, 0.0001)
    var"alpha.Age" ~ dnorm(0.0, 0.0001)
    var"alpha.V4" ~ dnorm(0.0, 0.0001)
    var"tau.b1" ~ dgamma(0.001, 0.001)
    var"sigma.b1" = 1.0 / sqrt(var"tau.b1")
    var"tau.b" ~ dgamma(0.001, 0.001)
    var"sigma.b" = 1.0 / sqrt(var"tau.b")
    alpha0 = ((((a0 - var"alpha.Base" * var"log.Base4.bar") - var"alpha.Trt" * var"Trt.bar") - var"alpha.BT" * var"BT.bar") - var"alpha.Age" * var"log.Age.bar") - var"alpha.V4" * var"V4.bar"
end

Data

The data set is too large to display comfortably here, so we load it from the copy that ships with JuliaBUGS. It contains N = 59 patients and T = 4 visits, the 59 × 4 matrix y of seizure counts, the treatment indicator Trt (0 = placebo, 1 = active treatment), the baseline seizure count Base, each patient's Age in years, and V4, an indicator that equals 1 only at the fourth visit.

data  = JuliaBUGS.BUGSExamples.VOLUME_1.epil.data
inits = JuliaBUGS.BUGSExamples.VOLUME_1.epil.inits
model = epil(data, inits)

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

  Model parameters:
    tau.b
    alpha.Base
    b[59, 4], b[59, 3], b[59, 2], b[59, 1], b[58, 4], b[58, 3], b[58, 2], b[58, 1], b[57, 4], b[57, 3], b[57, 2], b[57, 1], b[56, 4], b[56, 3], b[56, 2], b[56, 1], b[55, 4], b[55, 3], b[55, 2], b[55, 1], b[54, 4], b[54, 3], b[54, 2], b[54, 1], b[53, 4], b[53, 3], b[53, 2], b[53, 1], b[52, 4], b[52, 3], b[52, 2], b[52, 1], b[51, 4], b[51, 3], b[51, 2], b[51, 1], b[50, 4], b[50, 3], b[50, 2], b[50, 1], b[49, 4], b[49, 3], b[49, 2], b[49, 1], b[48, 4], b[48, 3], b[48, 2], b[48, 1], b[47, 4], b[47, 3], b[47, 2], b[47, 1], b[46, 4], b[46, 3], b[46, 2], b[46, 1], b[45, 4], b[45, 3], b[45, 2], b[45, 1], b[44, 4], b[44, 3], b[44, 2], b[44, 1], b[43, 4], b[43, 3], b[43, 2], b[43, 1], b[42, 4], b[42, 3], b[42, 2], b[42, 1], b[41, 4], b[41, 3], b[41, 2], b[41, 1], b[40, 4], b[40, 3], b[40, 2], b[40, 1], b[39, 4], b[39, 3], b[39, 2], b[39, 1], b[38, 4], b[38, 3], b[38, 2], b[38, 1], b[37, 4], b[37, 3], b[37, 2], b[37, 1], b[36, 4], b[36, 3], b[36, 2], b[36, 1], b[35, 4], b[35, 3], b[35, 2], b[35, 1], b[34, 4], b[34, 3], b[34, 2], b[34, 1], b[33, 4], b[33, 3], b[33, 2], b[33, 1], b[32, 4], b[32, 3], b[32, 2], b[32, 1], b[31, 4], b[31, 3], b[31, 2], b[31, 1], b[30, 4], b[30, 3], b[30, 2], b[30, 1], b[29, 4], b[29, 3], b[29, 2], b[29, 1], b[28, 4], b[28, 3], b[28, 2], b[28, 1], b[27, 4], b[27, 3], b[27, 2], b[27, 1], b[26, 4], b[26, 3], b[26, 2], b[26, 1], b[25, 4], b[25, 3], b[25, 2], b[25, 1], b[24, 4], b[24, 3], b[24, 2], b[24, 1], b[23, 4], b[23, 3], b[23, 2], b[23, 1], b[22, 4], b[22, 3], b[22, 2], b[22, 1], b[21, 4], b[21, 3], b[21, 2], b[21, 1], b[20, 4], b[20, 3], b[20, 2], b[20, 1], b[19, 4], b[19, 3], b[19, 2], b[19, 1], b[18, 4], b[18, 3], b[18, 2], b[18, 1], b[17, 4], b[17, 3], b[17, 2], b[17, 1], b[16, 4], b[16, 3], b[16, 2], b[16, 1], b[15, 4], b[15, 3], b[15, 2], b[15, 1], b[14, 4], b[14, 3], b[14, 2], b[14, 1], b[13, 4], b[13, 3], b[13, 2], b[13, 1], b[12, 4], b[12, 3], b[12, 2], b[12, 1], b[11, 4], b[11, 3], b[11, 2], b[11, 1], b[10, 4], b[10, 3], b[10, 2], b[10, 1], b[9, 4], b[9, 3], b[9, 2], b[9, 1], b[8, 4], b[8, 3], b[8, 2], b[8, 1], b[7, 4], b[7, 3], b[7, 2], b[7, 1], b[6, 4], b[6, 3], b[6, 2], b[6, 1], b[5, 4], b[5, 3], b[5, 2], b[5, 1], b[4, 4], b[4, 3], b[4, 2], b[4, 1], b[3, 4], b[3, 3], b[3, 2], b[3, 1], b[2, 4], b[2, 3], b[2, 2], b[2, 1], b[1, 4], b[1, 3], b[1, 2], b[1, 1]
    alpha.BT
    alpha.V4
    alpha.Age
    a0
    b1[59], b1[58], b1[57], b1[56], b1[55], b1[54], b1[53], b1[52], b1[51], b1[50], b1[49], b1[48], b1[47], b1[46], b1[45], b1[44], b1[43], b1[42], b1[41], b1[40], b1[39], b1[38], b1[37], b1[36], b1[35], b1[34], b1[33], b1[32], b1[31], b1[30], b1[29], b1[28], b1[27], b1[26], b1[25], b1[24], b1[23], b1[22], b1[21], b1[20], b1[19], b1[18], b1[17], b1[16], b1[15], b1[14], b1[13], b1[12], b1[11], b1[10], b1[9], b1[8], b1[7], b1[6], b1[5], b1[4], b1[3], b1[2], b1[1]
    tau.b1
    alpha.Trt

  Variable sizes and types:
    alpha.Base: type = Int64
    b: size = (59, 4), type = Matrix{Float64}
    alpha.BT: type = Int64
    Trt.bar: type = Float64
    V4: size = (4,), type = Vector{Int64}
    sigma.b1: type = Float64
    log.Base4: size = (59,), type = Vector{Float64}
    alpha.Trt: type = Int64
    tau.b: type = Int64
    V4.bar: type = Float64
    Age: size = (59,), type = Vector{Int64}
    alpha.V4: type = Int64
    log.Base4.bar: type = Float64
    Base: size = (59,), type = Vector{Int64}
    BT: size = (59,), type = Vector{Float64}
    b1: size = (59,), type = Vector{Float64}
    alpha0: type = Float64
    Trt: size = (59,), type = Vector{Int64}
    N: type = Int64
    BT.bar: type = Float64
    sigma.b: type = Float64
    mu: size = (59, 4), type = Matrix{Float64}
    T: type = Int64
    y: size = (59, 4), type = Matrix{Int64}
    log.Age.bar: type = Float64
    log.Age: size = (59,), type = Vector{Float64}
    alpha.Age: type = Int64
    a0: type = Int64
    tau.b1: type = Int64

We pass the example's published initial values as the second argument to epil rather than calling epil(data) alone. This is a log-linear Poisson model with vague priors, so values drawn at random from the priors can produce an invalid (non-positive or overflowing) rate; starting from sensible values keeps construction and sampling stable.

All of the classic Volume 1 examples ship with the package in JuliaBUGS.BUGSExamples, each providing the model definition, data, initial values, and reference results.

Sampling

To draw posterior samples, we build the model with gradient support and run the NUTS sampler from AdvancedHMC:

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

model = epil(data, inits; 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)

BUGS-style initial values for this example are available as JuliaBUGS.BUGSExamples.VOLUME_1.epil.inits and can be applied with initialize!(model, inits).

Results

JuliaBUGS does not ship reference results for this example. For published estimates, see the OpenBUGS documentation page, which reports results alongside the approximate-likelihood fit of Breslow and Clayton (1993); a correctly converged chain's summarystats should agree with those values up to Monte Carlo error.

See also: the example gallery overview and the getting started tutorial.