[ Note to Martyn Plummer : this arises from a thread on the lme4
development list : I tried to use JAGS to solve a mixed-model related
problem, and the (non-)results seems fishy... I the took the liberty to
Cc you on this subject ]
Le jeudi 04 mars 2010 ? 12:06 -0600, Douglas Bates a ?crit :
On Thu, Mar 4, 2010 at 11:32 AM, Federico Calboli
<f.calboli at imperial.ac.uk> wrote:
On 4 Mar 2010, at 16:11, Douglas Bates wrote:
<cut>
which, now I see are more consistent with the results from lme4, not
lme4a. I misunderstood the message I received regarding the Stata
results.
OK, this might have all been a red herring. I'll look into it some more.
At one time on Saturday Night Live Gilda Radner would play a citizen
commentator on weekend news update who was confused about the issue
and, after a long harangue, would end up realizing she got it wrong.
She ended by saying "Never mind". That might be the case here too.
Just out of curiosity, how do we know stata and lme4 are correct and lme4a is not? why is stata the benchmark? Please note I am only being facetious.
In cases like this we usually go with majority rule. I'm still
verifying and validating.
I'm currently trying to fit this in BUGS (JAGS, actually), and run in
strange difficulties : OpenBUGS (3.0.7) starts an endless loop while
adapting (first update after initialization), WinBUGS (1.4.3 with patch
in Wine) stops immediately with an unhelpful "Rejection 1" message. JAGS
(1.0.3) standalone seems to load model, load data, initializing well.
However, it stops at the adaptation step, stating "Error in node
alpha[<variable...>]. Current value is inconsistent with data". JAGS
called from R (through Martyn Plummer's rjags package) *crashes R*
(almost never seen that before...). I get this in ESS :
terminate called after throwing an instance of 'std::out_of_range'
what(): Range: upper < lower bound in constructor
Process R abandon at Thu Mar 4 20:37:46 2010
The model, in which I find no "obvious" fault, is as follows :
model
{
for (k in 1:nobs) {
incidence[k] ~ dbin(p[k], size[k])
logit(p[k]) <- alpha[herd[k]] + beta[period[k]] + gamma
}
gamma ~ dnorm(0.00000E+00, tau.gamma)
tau.gamma <- pow(sigma.gamma, -2)
sigma.gamma ~ dunif(0.00000E+00, 100)
for (i in 1:nherd) {
alpha[i] ~ dnorm(0.00000E+00, tau.alpha)
}
tau.alpha <- pow(sigma.alpha, -2)
sigma.alpha ~ dunif(0.00000E+00, 100)
for (j in 1:nperiod) {
beta[j] ~ dnorm(0.00000E+00, 1.00000E-06)
}
}
I tried to pass "reasonable" initial values through model.jags (rjags) :
inits=function() {
list(mu.alpha=rnorm(nlevels(cbpp$herd),0,1),
mu.beta=rnorm(nlevels(cbpp
$period),0,1),
gamma=rnorm(0,1))
... to no avail ! I even took the usual "unusual step" of looking at the
raw data and graphing them, with nothing suspicious jumping out of the
graphs...
I'm Cc'ing Martyn Plummer on this one, and am curious of his thoughts on
the subject...
HTH,
Emmanuel Charpentier