GitHub

Oxford: Smooth Fit to Log-Odds Ratios

This example comes from a classic case-control study of childhood cancer and maternal exposure to X-rays during pregnancy, analysed by Breslow and Clayton (1993). The data are arranged as 120 small 2-by-2 tables, one per stratum. Each stratum cross-classifies cases (children who died of cancer) and matched controls by whether the mother received a prenatal X-ray, and the strata are defined by the child's age group and birth-year cohort spanning the years 1944 to 1964. For stratum $i$, r1[i] of the n1[i] cases were exposed and r0[i] of the n0[i] controls were exposed, while year[i] records the (centred) birth year of that stratum.

The scientific question is whether the association between prenatal X-ray exposure and childhood cancer changed smoothly over the birth years covered by the study. The model is a Bayesian hierarchical binomial (logistic) regression: each stratum gets its own nuisance intercept mu[i] for the exposure odds among controls, and the log-odds ratio comparing cases with controls, logPsi[i], is described as a smooth quadratic function of birth year plus a normal random effect that absorbs residual between-stratum variation. It is one of the examples in Volume 1 of the classic BUGS examples; see also the OpenBUGS version of this example.

Model

\[\begin{aligned} r0_i &\sim \text{Binomial}(n0_i,\ p0_i), & \text{logit}(p0_i) &= \mu_i \\ r1_i &\sim \text{Binomial}(n1_i,\ p1_i), & \text{logit}(p1_i) &= \mu_i + \log\!\Psi_i \\ \log\!\Psi_i &= \alpha + \beta_1\, \text{year}_i + \beta_2\, (\text{year}_i^2 - 22) + b_i \\ b_i &\sim \text{Normal}(0, \tau) \\ \mu_i &\sim \text{Normal}(0, 10^{-6}) \end{aligned}\]

The population parameters $\alpha$, $\beta_1$, $\beta_2$ are given "noninformative" normal priors, $\tau$ a Gamma(0.001, 0.001) prior, and $\sigma = 1/\sqrt{\tau}$ is the random-effect standard deviation. Here $\tau$ denotes the precision (inverse variance) of a normal distribution, following the BUGS convention.

using JuliaBUGS

oxford = @bugs begin
    for i in 1:K
        r0[i] ~ dbin(p0[i], n0[i])
        r1[i] ~ dbin(p1[i], n1[i])
        p0[i] = logistic(mu[i])
        p1[i] = logistic(mu[i] + logPsi[i])
        logPsi[i] = alpha + beta1 * year[i] + beta2 * (year[i] * year[i] - 22) + b[i]
        b[i] ~ dnorm(0, tau)
        mu[i] ~ dnorm(0.0, 1.0e-6)
    end
    alpha ~ dnorm(0.0, 1.0e-6)
    beta1 ~ dnorm(0.0, 1.0e-6)
    beta2 ~ dnorm(0.0, 1.0e-6)
    tau ~ dgamma(0.001, 0.001)
    sigma = 1 / sqrt(tau)
end
BUGSModelDef:
begin
    for i = 1:K
        r0[i] ~ dbin(p0[i], n0[i])
        r1[i] ~ dbin(p1[i], n1[i])
        p0[i] = logistic(mu[i])
        p1[i] = logistic(mu[i] + logPsi[i])
        logPsi[i] = alpha + beta1 * year[i] + beta2 * (year[i] * year[i] - 22) + b[i]
        b[i] ~ dnorm(0, tau)
        mu[i] ~ dnorm(0.0, 1.0e-6)
    end
    alpha ~ dnorm(0.0, 1.0e-6)
    beta1 ~ dnorm(0.0, 1.0e-6)
    beta2 ~ 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 with the number of strata K, the centred birth year year of each stratum, and the exposure counts among cases (r1 out of n1) and among controls (r0 out of n0).

data = (
    r1 = [3, 5, 2, 7, 7, 2, 5, 3, 5, 11, 6, 6, 11, 4, 4, 2, 8, 8, 6, 5, 15, 4, 9, 9, 4,
        12, 8, 8, 6, 8, 12, 4, 7, 16, 12, 9, 4, 7, 8, 11, 5, 12, 8, 17, 9, 3, 2, 7, 6,
        5, 11, 14, 13, 8, 6, 4, 8, 4, 8, 7, 15, 15, 9, 9, 5, 6, 3, 9, 12, 14, 16, 17,
        8, 8, 9, 5, 9, 11, 6, 14, 21, 16, 6, 9, 8, 9, 8, 4, 11, 11, 6, 9, 4, 4, 9, 9,
        10, 14, 6, 3, 4, 6, 10, 4, 3, 3, 10, 4, 10, 5, 4, 3, 13, 1, 7, 5, 7, 6, 3, 7],
    n1 = [28, 21, 32, 35, 35, 38, 30, 43, 49, 53, 31, 35, 46, 53, 61, 40, 29, 44, 52, 55,
        61, 31, 48, 44, 42, 53, 56, 71, 43, 43, 43, 40, 44, 70, 75, 71, 37, 31, 42, 46,
        47, 55, 63, 91, 43, 39, 35, 32, 53, 49, 75, 64, 69, 64, 49, 29, 40, 27, 48, 43,
        61, 77, 55, 60, 46, 28, 33, 32, 46, 57, 56, 78, 58, 52, 31, 28, 46, 42, 45, 63,
        71, 69, 43, 50, 31, 34, 54, 46, 58, 62, 52, 41, 34, 52, 63, 59, 88, 62, 47, 53,
        57, 74, 68, 61, 45, 45, 62, 73, 53, 39, 45, 51, 55, 41, 53, 51, 42, 46, 54, 32],
    r0 = [0, 2, 2, 1, 2, 0, 1, 1, 1, 2, 4, 4, 2, 1, 7, 4, 3, 5, 3, 2, 4, 1, 4, 5, 2,
        7, 5, 8, 2, 3, 5, 4, 1, 6, 5, 11, 5, 2, 5, 8, 5, 6, 6, 10, 7, 5, 5, 2, 8,
        1, 13, 9, 11, 9, 4, 4, 8, 6, 8, 6, 8, 14, 6, 5, 5, 2, 4, 2, 9, 5, 6, 7,
        5, 10, 3, 2, 1, 7, 9, 13, 9, 11, 4, 8, 2, 3, 7, 4, 7, 5, 6, 6, 5, 6, 9, 7,
        7, 7, 4, 2, 3, 4, 10, 3, 4, 2, 10, 5, 4, 5, 4, 6, 5, 3, 2, 2, 4, 6, 4, 1],
    n0 = [28, 21, 32, 35, 35, 38, 30, 43, 49, 53, 31, 35, 46, 53, 61, 40, 29, 44, 52, 55,
        61, 31, 48, 44, 42, 53, 56, 71, 43, 43, 43, 40, 44, 70, 75, 71, 37, 31, 42, 46,
        47, 55, 63, 91, 43, 39, 35, 32, 53, 49, 75, 64, 69, 64, 49, 29, 40, 27, 48, 43,
        61, 77, 55, 60, 46, 28, 33, 32, 46, 57, 56, 78, 58, 52, 31, 28, 46, 42, 45, 63,
        71, 69, 43, 50, 31, 34, 54, 46, 58, 62, 52, 41, 34, 52, 63, 59, 88, 62, 47, 53,
        57, 74, 68, 61, 45, 45, 62, 73, 53, 39, 45, 51, 55, 41, 53, 51, 42, 46, 54, 32],
    year = [
        -10, -9, -9, -8, -8, -8, -7, -7, -7, -7, -6, -6, -6, -6, -6, -5, -5, -5, -5, -5, -5,
        -4, -4, -4, -4, -4, -4, -4, -3, -3, -3, -3, -3, -3, -3, -3, -2, -2, -2, -2, -2, -2,
        -2, -2, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4,
        4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 9, 9, 10],
    K = 120
)

model = oxford(data)
BUGSModel (parameters are in transformed (unconstrained) space, with dimension 244):

  Model parameters:
    alpha
    beta2
    b[120], b[119], b[118], b[117], b[116], b[115], b[114], b[113], b[112], b[111], b[110], b[109], b[108], b[107], b[106], b[105], b[104], b[103], b[102], b[101], b[100], b[99], b[98], b[97], b[96], b[95], b[94], b[93], b[92], b[91], b[90], b[89], b[88], b[87], b[86], b[85], b[84], b[83], b[82], b[81], b[80], b[79], b[78], b[77], b[76], b[75], b[74], b[73], b[72], b[71], b[70], b[69], b[68], b[67], b[66], b[65], b[64], b[63], b[62], b[61], b[60], b[59], b[58], b[57], b[56], b[55], b[54], b[53], b[52], b[51], b[50], b[49], b[48], b[47], b[46], b[45], b[44], b[43], b[42], b[41], b[40], b[39], b[38], b[37], b[36], b[35], b[34], b[33], b[32], b[31], b[30], b[29], b[28], b[27], b[26], b[25], b[24], b[23], b[22], 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]
    mu[120], mu[119], mu[118], mu[117], mu[116], mu[115], mu[114], mu[113], mu[112], mu[111], mu[110], mu[109], mu[108], mu[107], mu[106], mu[105], mu[104], mu[103], mu[102], mu[101], mu[100], mu[99], mu[98], mu[97], mu[96], mu[95], mu[94], mu[93], mu[92], mu[91], mu[90], mu[89], mu[88], mu[87], mu[86], mu[85], mu[84], mu[83], mu[82], mu[81], mu[80], mu[79], mu[78], mu[77], mu[76], mu[75], mu[74], mu[73], mu[72], mu[71], mu[70], mu[69], mu[68], mu[67], mu[66], mu[65], mu[64], mu[63], mu[62], mu[61], mu[60], mu[59], mu[58], mu[57], mu[56], mu[55], mu[54], mu[53], mu[52], mu[51], mu[50], mu[49], mu[48], mu[47], mu[46], mu[45], mu[44], mu[43], mu[42], mu[41], mu[40], mu[39], mu[38], mu[37], mu[36], mu[35], mu[34], mu[33], mu[32], mu[31], mu[30], mu[29], mu[28], mu[27], mu[26], mu[25], mu[24], mu[23], mu[22], mu[21], mu[20], mu[19], mu[18], mu[17], mu[16], mu[15], mu[14], mu[13], mu[12], mu[11], mu[10], mu[9], mu[8], mu[7], mu[6], mu[5], mu[4], mu[3], mu[2], mu[1]
    tau
    beta1

  Variable sizes and types:
    alpha: type = Float64
    b: size = (120,), type = Vector{Float64}
    n1: size = (120,), type = Vector{Int64}
    p0: size = (120,), type = Vector{Float64}
    p1: size = (120,), type = Vector{Float64}
    n0: size = (120,), type = Vector{Int64}
    year: size = (120,), type = Vector{Int64}
    r1: size = (120,), type = Vector{Int64}
    K: type = Int64
    sigma: type = Float64
    mu: size = (120,), type = Vector{Float64}
    logPsi: size = (120,), type = Vector{Float64}
    tau: type = Float64
    beta1: type = Float64
    r0: size = (120,), type = Vector{Int64}
    beta2: type = Float64

All of the classic examples ship with the package in JuliaBUGS.BUGSExamples, so the model definition, data, initial values, and reference results are also available directly as JuliaBUGS.BUGSExamples.VOLUME_1.oxford.

Sampling

We draw posterior samples with the NUTS sampler from AdvancedHMC, rebuilding the model with gradient support first.

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

model = oxford(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)

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

Results

This transcription does not bundle a numeric reference table (JuliaBUGS.BUGSExamples.VOLUME_1.oxford.reference_results is nothing), so there is nothing to reproduce verbatim here. For published posterior summaries of $\alpha$, $\beta_1$, $\beta_2$, and $\sigma$, consult the OpenBUGS Oxford example and the original analysis by Breslow and Clayton (1993). A correctly converged chain's summarystats should agree with those published values up to Monte Carlo error.

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