[R-sig-ME] Cross-validated likelihood, cont.
d@r|zopou|o@ @end|ng |rom er@@mu@mc@n|
Mon Apr 29 06:25:01 CEST 2019
If you want you could give a try to the GLMMadaptive package that implements the adaptive Gaussian quadrature for a vector of random effects (e.g., intercepts and slopes as in your case), and from which you get the log-likelihood in two steps, e.g.,
# the log-likelihood at the initial values
fm <- mixed_model(cbind(Dead, Alive) ~ (Trt + 0) / Dose, random = ~ Dose | Rep,
data = Ts, family = binomial(link = �cloglog�), iter_EM = 0, iter_qN_outer = 0)
# the log-likelihood at user-specified values
gm <- mixed_model(cbind(Dead, Alive) ~ (Trt + 0) / Dose, random = ~ Dose | Rep,
data = Ts, family = binomial(link = �cloglog�), iter_EM = 0, iter_qN_outer = 0,
initial_values = list(beta = <put_your_fixed_effects_here>, D = <put_the_RE_cov_matrix_here>))
- - - - - -
Professor of Biostatistics
Erasmus University Medical Center
From: Rolf Turner <r.turner using auckland.ac.nz<mailto:r.turner using auckland.ac.nz>>
Date: Monday, 29 Apr 2019, 2:41 AM
To: Ben Bolker <bbolker using gmail.com<mailto:bbolker using gmail.com>>
Cc: r-sig-mixed-models using r-project.org <r-sig-mixed-models using r-project.org<mailto:r-sig-mixed-models using r-project.org>>
Subject: Re: [R-sig-ME] Cross-validated likelihood, cont.
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),
coefs <- unlist(getME(f.trn,"theta"))
newdev <- update(f.trn, data=VS, devFunOnly=TRUE,
check <- update(f.trn, data=TS, devFunOnly=TRUE,
which are way out of line with the nAGQ=0 results:
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
Or am I, as is so often the case ( :-) ) out of luck?
> 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:
>>> [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
>>> 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),
>>> ## 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
>>> newdev <- update(f.trn, data=VS, devFunOnly=TRUE,
R-sig-mixed-models using r-project.org mailing list
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models