[R-sig-ME] Cross-validated likelihood, cont.
D. Rizopoulos
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.,
library(GLMMadaptive)
# 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)
logLik(fm)
# 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>))
logLik(gm)
Best,
Dimitris
- - - - - -
Dimitris Rizopoulos
Professor of Biostatistics
Erasmus University Medical Center
The Netherlands
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),
data=TS,family=binomial(link="cloglog"),nAGQ=1)
coefs <- unlist(getME(f.trn,"theta"))
newdev <- update(f.trn, data=VS, devFunOnly=TRUE,
control=glmerControl(check.nobs.vs.nRE="ignore"))
print(newdev(coefs))
check <- update(f.trn, data=TS, devFunOnly=TRUE,
control=glmerControl(check.nobs.vs.nRE="ignore"))
print(check(coefs))
I get:
[1] 390.7528
and
[1] 2194.427
which are way out of line with the nAGQ=0 results:
[1] 42.75071
and
[1] 245.8587
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
nAGQ=0 fits.)
Or am I, as is so often the case ( :-) ) out of luck?
cheers,
Rolf
>
> 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)
_______________________________________________
R-sig-mixed-models using r-project.org mailing list
https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-mixed-models&data=02%7C01%7Cd.rizopoulos%40erasmusmc.nl%7C2809ad114d6c4906bfec08d6cc330a1e%7C526638ba6af34b0fa532a1a511f4ac80%7C0%7C0%7C636920916920782143&sdata=nds0cDpFmm2PAcrP1g4NlhabC9u85taAdUe2u213UGo%3D&reserved=0
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list