[R-sig-ME] Cross-validated likelihood, cont.
Rolf Turner
r@turner @end|ng |rom @uck|@nd@@c@nz
Mon Apr 29 01:03:18 CEST 2019
Please see in-line below.
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).
Some rudimentary experimentation that I've done appears to indicate that
the *number* returned is equal to -2*(the log likelihood) --- presumably
the deviance. Is this correct?
If I want to get the log likelihood (e.g. for comparison with other
procedures) is it "safe" to assume that the log likelihood is equal to
-0.5*(the number returned)?
Thanks.
cheers,
Rolf
>
> 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