LeukFr: Cox Regression with Random Effects

This example comes from Volume 1 of the classic BUGS examples, and is also documented on the OpenBUGS Leukfr page. The data are the well-known leukaemia remission times of Freireich et al. (1963): 42 patients arranged in 21 matched pairs, where within each pair one patient received the drug 6-mercaptopurine (6-MP) and the other received placebo. For each patient we observe a remission (survival) time obs.t and an indicator fail of whether the event was observed or the observation was censored, together with a treatment covariate Z (coded as ±0.5) and the pairing pair.

The question is whether 6-MP prolongs remission, while accounting for the fact that patients were matched in pairs. The model is a Cox proportional-hazards regression fitted through the counting-process (Poisson) representation used throughout the BUGS survival examples. It extends the plain Leuk model by adding a normally distributed frailty term b[pair[i]], a random effect shared by the two patients in each matched pair, so that within-pair correlation is modelled explicitly. This is a random-effects (frailty) survival model.

Model

At each observed failure time t[j] the counting-process increment dN[i, j] for patient i is treated as Poisson with an intensity that combines the treatment effect, the pair-specific frailty, and a baseline hazard increment dL0[j]:

\[\begin{aligned} dN_{ij} &\sim \text{Poisson}(I_{ij}) \\ I_{ij} &= Y_{ij}\,\exp(\beta\, Z_i + b_{\text{pair}(i)})\,dL0_j \\ b_k &\sim \text{Normal}(0, \tau) \\ dL0_j &\sim \text{Gamma}(\mu_j, c) \\ \beta &\sim \text{Normal}(0, 10^{-6}) \\ \tau &\sim \text{Gamma}(0.001, 0.001) \end{aligned}\]

Here Y[i, j] is the risk-set indicator (1 if patient i is still at risk at time t[j]), built from the observed times using step. Some variables in the original program have R-style dotted names such as obs.t, dL0.star, S.treat, and S.placebo; these are written with Julia's var"..." syntax so the names carry over from the original BUGS program unchanged.

using JuliaBUGS

leukfr = @bugs begin # Set up data
    for i in 1:N
        for j in 1:T
            # risk set = 1 if obs.t >= t
            Y[i, j] = step(var"obs.t"[i] - t[j] + eps)

            # counting process jump = 1 if obs.t in [ t[j], t[j+1] )
            # i.e. if t[j] <= obs.t < t[j+1]
            dN[i, j] = Y[i, j] * step(t[j + 1] - var"obs.t"[i] - eps) * fail[i]
        end
    end

    # Model
    for j in 1:T
        for i in 1:N
            dN[i, j] ~ dpois(Idt[i, j])
            Idt[i, j] = Y[i, j] * exp(beta * Z[i] + b[pair[i]]) * dL0[j]
        end
        dL0[j] ~ dgamma(mu[j], c)
        mu[j] = var"dL0.star"[j] * c # prior mean hazard

        # Survivor function = exp(-Integral{l0(u)du})^exp(beta * z)
        var"S.treat"[j] = pow(exp(-sum(dL0[1:j])), exp(beta * -0.5))
        var"S.placebo"[j] = pow(exp(-sum(dL0[1:j])), exp(beta * 0.5))
    end
    for k in 1:Npairs
        b[k] ~ dnorm(0.0, tau)
    end
    tau ~ dgamma(0.001, 0.001)
    sigma = sqrt(1 / tau)
    c = 0.001
    r = 0.1
    for j in 1:T
        var"dL0.star"[j] = r * (t[j + 1] - t[j])
    end
    beta ~ dnorm(0.0, 0.000001)
end

BUGSModelDef:
begin
    for i = 1:N
        for j = 1:T
            Y[i, j] = _step((var"obs.t"[i] - t[j]) + eps)
            dN[i, j] = Y[i, j] * _step((t[j + 1] - var"obs.t"[i]) - eps) * fail[i]
        end
    end
    for j = 1:T
        for i = 1:N
            dN[i, j] ~ dpois(Idt[i, j])
            Idt[i, j] = Y[i, j] * exp(beta * Z[i] + b[pair[i]]) * dL0[j]
        end
        dL0[j] ~ dgamma(mu[j], c)
        mu[j] = var"dL0.star"[j] * c
        var"S.treat"[j] = pow(exp(-(sum(dL0[1:j]))), exp(beta * -0.5))
        var"S.placebo"[j] = pow(exp(-(sum(dL0[1:j]))), exp(beta * 0.5))
    end
    for k = 1:Npairs
        b[k] ~ dnorm(0.0, tau)
    end
    tau ~ dgamma(0.001, 0.001)
    sigma = sqrt(1 / tau)
    c = 0.001
    r = 0.1
    for j = 1:T
        var"dL0.star"[j] = r * (t[j + 1] - t[j])
    end
    beta ~ dnorm(0.0, 1.0e-6)
end

Data

The data give the study dimensions (N = 42 patients, T = 17 distinct failure times, Npairs = 21 matched pairs), the grid of failure times t, the observed remission times obs.t, the failure/censoring indicators fail, the pair labels pair, and the treatment covariate Z. We supply them as a NamedTuple and construct the model:

data = (
    N = 42,
    T = 17,
    eps = 0.00001,
    Npairs = 21,
    t = [1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 13, 15, 16, 17, 22, 23, 35],
    var"obs.t" = [
        1, 1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 8, 8, 11, 11, 12, 12, 15, 17, 22, 23, 6,
        6, 6, 6, 7, 9, 10, 10, 11, 13, 16, 17, 19, 20, 22, 23, 25, 32, 32, 34, 35
    ],
    pair = [
        1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
        19, 18, 8, 1, 20, 6, 2, 10, 3, 14, 4, 11, 7, 9, 12, 16, 17, 5, 13, 15, 21
    ],
    fail = [
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0
    ],
    Z = [
        0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
        0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
        0.5, 0.5, 0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
        -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
        -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5
    ]
)

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

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

  Model parameters:
    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]
    beta
    dL0[17], dL0[16], dL0[15], dL0[14], dL0[13], dL0[12], dL0[11], dL0[10], dL0[9], dL0[8], dL0[7], dL0[6], dL0[5], dL0[4], dL0[3], dL0[2], dL0[1]
    tau

  Variable sizes and types:
    b: size = (21,), type = Vector{Float64}
    dL0.star: size = (17,), type = Vector{Float64}
    c: type = Float64
    eps: type = Float64
    dN: size = (42, 17), type = Matrix{Float64}
    obs.t: size = (42,), type = Vector{Int64}
    S.treat: size = (17,), type = Vector{Float64}
    Npairs: type = Int64
    S.placebo: size = (17,), type = Vector{Float64}
    dL0: size = (17,), type = Vector{Float64}
    sigma: type = Float64
    Idt: size = (42, 17), type = Matrix{Float64}
    pair: size = (42,), type = Vector{Int64}
    Z: size = (42,), type = Vector{Float64}
    N: type = Int64
    fail: size = (42,), type = Vector{Int64}
    mu: size = (17,), type = Vector{Float64}
    Y: size = (42, 17), type = Matrix{Float64}
    tau: type = Int64
    r: type = Float64
    T: type = Int64
    beta: type = Float64
    t: size = (18,), type = Vector{Int64}

We pass the example's published initial values as the second argument to leukfr rather than calling leukfr(data) alone: this counting-process (Poisson) survival model has vague priors, so values drawn at random from the priors can give an invalid rate, and a sensible starting point keeps construction and sampling stable.

All of the classic examples ship with the package under JuliaBUGS.BUGSExamples, each bundling the model definition, data, initial values, and reference results, so you can also load this one directly as JuliaBUGS.BUGSExamples.VOLUME_1.leukfr.

Sampling

Construct the model with a gradient backend and draw posterior samples with the No-U-Turn sampler:

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

model = leukfr(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.leukfr.inits and can be applied with initialize!(model, inits).

Results

The source file ships reference_results = nothing for this example, so there is no tabulated reference posterior summary bundled with the package to reproduce here. For published numerical summaries of the treatment effect beta and the frailty standard deviation sigma (obtained after a 1,000-iteration burn-in followed by 10,000 further updates), see the OpenBUGS Leukfr page. A correctly converged chain's summarystats output should match those published values up to Monte Carlo error.

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