[R-sig-ME] Cross-validated likelihood, cont.

[Ben Bolker]

No, nAGQ=0 is probably as simple as you can get.

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)
>