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
where are the number of seeds, out of , that germinate on plate ; and and 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