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

Rolf Turner r@turner @end|ng |rom @uck|@nd@@c@nz
Fri Apr 26 06:41:00 CEST 2019


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