[R-sig-ME] Cross-validated likelihood, cont.
bbo|ker @end|ng |rom gm@||@com
Fri Apr 26 13:01:04 CEST 2019
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.
> 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