[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