[R-sig-ME] Cross-validated likelihood, cont.
Rolf Turner
r@turner @end|ng |rom @uck|@nd@@c@nz
Mon Apr 29 01:41:01 CEST 2019
On 26/04/19 11:01 PM, Ben Bolker wrote:
> No, nAGQ=0 is probably as simple as you can get.
Hmm. Indeed. I tried things just now with nAGQ=1, and everything seems
to go to custard. If I do:
f.trn <- glmer(cbind(Dead,Alive) ~ (Trt+0)/Dose + (Dose | Rep),
data=TS,family=binomial(link="cloglog"),nAGQ=1)
coefs <- unlist(getME(f.trn,"theta"))
newdev <- update(f.trn, data=VS, devFunOnly=TRUE,
control=glmerControl(check.nobs.vs.nRE="ignore"))
print(newdev(coefs))
check <- update(f.trn, data=TS, devFunOnly=TRUE,
control=glmerControl(check.nobs.vs.nRE="ignore"))
print(check(coefs))
I get:
[1] 390.7528
and
[1] 2194.427
which are way out of line with the nAGQ=0 results:
[1] 42.75071
and
[1] 245.8587
Note that when nAGQ is set equal to 0, logLik(f.trn) gives -122.9294
which is, at least approximately equal to -0.5*245.8587.
Note also that when nAGQ is set equal to 1, logLik(f.trn) gives
-122.8933 which is "reasonably" close to the nAGQ=0 result. Whereas
2194.427 is *not* by any stretch close to 245.8587.
Is there anything that can be done to make nAGQ=1 "behave"? Ideally I
would like to be able to apply "cross-validated" log likelihood to
models that are fitted any which-way. (And *not* have to restrict to
nAGQ=0 fits.)
Or am I, as is so often the case ( :-) ) out of luck?
cheers,
Rolf
>
> On Fri, Apr 26, 2019 at 12:41 AM Rolf Turner <r.turner using auckland.ac.nz> wrote:
>>
>>
>> Got it. Thanks very much. A bit messier than I'd hoped, but I guess
>> that's life ....
>>
>> I shall look at the lme4::modular example that you pointed to, and see
>> if it enlightens me at all.
>>
>> Would things have less of a tendency to go to custard if I set nAGQ=1?
>> I've been using nAGQ=0 since it's faster and less likely to throw
>> errors, and appears to give essentially the same results in the
>> applications in which I am interested.
>>
>> cheers,
>>
>> Rolf
>>
>> On 26/04/19 1:32 PM, Ben Bolker wrote:
>>> Rolf,
>>>
>>> [cc'ing r-sig-mixed-models; apologies for breaking the thread]
>>>
>>> Two issues:
>>>
>>> - I had to override some checking that glmer does (here I think the
>>> problem is that the sample was sufficiently unbalanced that glmer
>>> figured the random effects would be unidentifiable)
>>> - there's some trickiness about whether the deviance function is
>>> stage-1/nAGQ=0 -- in which case, as Doug pointed out, the fixed-effect
>>> [beta] parameters are profiled out, and the parameter vector should only
>>> include theta (var-cov/Cholesky params), not theta+beta -- or if it's
>>> stage-2/nAGQ>0 -- in which case the parameter vector should combine
>>> theta and beta.
>>>
>>> What's below seems to work (at least, it returns a number).
>>>
>>> If this seems mysterious, the second example in ?lme4::modular *might*
>>> help ...
>>>
>>>
>>> ====
>>> #
>>> # Script demo.txt
>>> #
>>>
>>> library(lme4)
>>> set.seed(42)
>>> X <- dget("X.txt")
>>> ind.trn <- sample(1:124,100)
>>> ind.val <- setdiff(1:124,ind.trn)
>>> TS <- X[ind.trn,]
>>> VS <- X[ind.val,]
>>> f.trn <- glmer(cbind(Dead,Alive) ~ (Trt+0)/Dose + (Dose | Rep),
>>> data=TS,family=binomial(link="cloglog"),nAGQ=0)
>>> ## coefs <- unlist(getME(f.trn,c("theta","beta")))
>>> ##
>>> coefs <- unlist(getME(f.trn,"theta"))
>>> ## Error: number of observations (=24) < number of random effects (=30)
>>> for term (Dose | Rep); the random-effects parameters are probably
>>> unidentifiable
>>> newdev <- update(f.trn, data=VS, devFunOnly=TRUE,
>>> control=glmerControl(check.nobs.vs.nRE="ignore"))
>>> newdev(coefs)
More information about the R-sig-mixed-models
mailing list