Seeds: Random Effect Logistic Regression

An example from OpenBUGS [43], Crowder [19], and Breslow and Clayton [10] concerning the proportion of seeds that germinated on each of 21 plates arranged according to a 2 by 2 factorial layout by seed and type of root extract.

Model

Germinations are modelled as

r_i &\sim \text{Binomial}(n_i, p_i) \quad\quad i=1,\ldots,21 \\
\operatorname{logit}(p_i) &= \alpha_0 + \alpha_1 x_{1i} + \alpha_2 x_{2i} + \alpha_{12} x_{1i} x_{2i} + b_i \\
b_i &\sim \text{Normal}(0, \sigma) \\
\alpha_0, \alpha_1, \alpha_2, \alpha_{12} &\sim \text{Normal}(0, 1000) \\
\sigma^2 &\sim \text{InverseGamma}(0.001, 0.001),

where r_i are the number of seeds, out of n_i, that germinate on plate i; and x_{1i} and x_{2i} are the seed type and root extract.

Analysis Program

using Mamba

## Data
seeds = Dict{Symbol, Any}(
  :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]
)
seeds[:N] = length(seeds[:r])


## Model Specification
model = Model(

  r = Stochastic(1,
    (alpha0, alpha1, x1, alpha2, x2, alpha12, b, n, N) ->
      UnivariateDistribution[
        begin
          p = invlogit(alpha0 + alpha1 * x1[i] + alpha2 * x2[i] +
                       alpha12 * x1[i] * x2[i] + b[i])
          Binomial(n[i], p)
        end
        for i in 1:N
      ],
    false
  ),

  b = Stochastic(1,
    s2 -> Normal(0, sqrt(s2)),
    false
  ),

  alpha0 = Stochastic(
    () -> Normal(0, 1000)
  ),

  alpha1 = Stochastic(
    () -> Normal(0, 1000)
  ),

  alpha2 = Stochastic(
    () -> Normal(0, 1000)
  ),

  alpha12 = Stochastic(
    () -> Normal(0, 1000)
  ),

  s2 = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  )

)


## Initial Values
inits = [
  Dict(:r => seeds[:r], :alpha0 => 0, :alpha1 => 0, :alpha2 => 0,
       :alpha12 => 0, :s2 => 0.01, :b => zeros(seeds[:N])),
  Dict(:r => seeds[:r], :alpha0 => 0, :alpha1 => 0, :alpha2 => 0,
       :alpha12 => 0, :s2 => 1, :b => zeros(seeds[:N]))
]


## Sampling Scheme
scheme = [AMM([:alpha0, :alpha1, :alpha2, :alpha12], 0.01 * eye(4)),
          AMWG(:b, 0.01),
          AMWG(:s2, 0.1)]
setsamplers!(model, scheme)


## MCMC Simulations
sim = mcmc(model, seeds, inits, 12500, burnin=2500, thin=2, chains=2)
describe(sim)

Results

Iterations = 2502:12500
Thinning interval = 2
Chains = 1,2
Samples per chain = 5000

Empirical Posterior Estimates:
            Mean         SD       Naive SE       MCSE        ESS
 alpha2  1.310728093 0.26053104 0.0026053104 0.0153996801 286.21707
 alpha1  0.088700176 0.26872879 0.0026872879 0.0128300598 438.70341
 alpha0 -0.556154341 0.17595432 0.0017595432 0.0101730837 299.15388
alpha12 -0.746440855 0.43006756 0.0043006756 0.0251658152 292.04607
     s2  0.085705306 0.09738014 0.0009738014 0.0080848189 145.07755

Quantiles:
             2.5%         25.0%        50.0%       75.0%       97.5%
 alpha2  0.8040593795  1.148881827  1.309947687  1.48076318  1.82815608
 alpha1 -0.4250164771 -0.093637900  0.094390643  0.26007581  0.62353933
 alpha0 -0.9149197759 -0.666632319 -0.551292851 -0.44262420 -0.22244775
alpha12 -1.5457041398 -1.027576522 -0.757250262 -0.49149187  0.17029702
     s2  0.0011739822  0.021117624  0.059376315  0.11140082  0.35234645