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

[Rolf Turner]

On 15/08/19 6:04 AM, D. Rizopoulos wrote:

> Dear Rolf,
> 
> The optimization behind mixed_model() is a hybrid that starts with EM
> and then switches to quasi-Newton. When you set the control arguments
> "iter_EM = 0" and "iter_qN_outer = 0" specify that neither EM nor
> quasi-Newton iterations are used. Hence, you only calculate the
> log-likelihood value (and the rest of the components) only for the
> initial values.
> 
> The reason why you received the error message is because in the list
> you provide to the initial_values argument of mixed_model() should be
> a named list. The name for the fixed effects coefficients is "betas"
> not "beta" - see the online help file for more info. I had also made
> this mistake in my untested code you quoted below. If you make this
> change, it should work (at least it does in my machine).

Yesssss! That worked!!!  Thank you.  And thank you for confirming that 
my understanding of what is going on was basically correct.

What a difference an "s" makes!  Especially when one hasn't a clue what 
one is doing. :-( Of course this could be held up as another instance of 
RTFM!  Now that I *look* at the help for mixed_model() I see that it 
clearly states that the components of "initial_values" are named 
*"betas"*, "D" and "phis".

I was feeling so all-at-sea that I never considered looking for such a 
simple solution to the problem.

Thanks again.

cheers,

Rolf

P.S.  (The following is mostly aimed at Ben Bolker.)

I thought I'd compare the log likelihood from mixed_model() with that 
from glmer().  So I did:

library(lme4)
g.trn   <- glmer(cbind(Dead,Alive) ~ (Trt+0)/Dose + (Dose | Rep),
                  data=TS,family=binomial(link="cloglog"),nAGQ=0)
coefs   <- unlist(getME(g.trn,c("theta","beta")))
newdev  <- update(g.trn, data=VS, devFunOnly=TRUE,
                   control=glmerControl(check.nobs.vs.nRE="ignore"))
ll <- -0.5*newdev(coefs)
print(ll)

This *worked* previously.  However after the "newdev <-" it now makes 
the comment:

> fixed-effect model matrix is rank deficient so dropping 1 column / coefficient

and then after "newdev(coefs) it throws an error:

> (15!=3)
> Error in pp$setTheta(theta) : theta size mismatch

I thought that the dropped coefficient might be the problem, so I 
replaced coefs by coefs[-15] but then the error simply changes to:

> (14!=3)
> Error in pp$setTheta(theta) : theta size mismatch

Indeed 15!=3 (no shit, Sherlock!) and likewise 14!=3.  But why *3*???

If I do

chk <- update(g.trn,data=VS)
coefs.chk <- unlist(getME(chk,c("theta","beta")))

I get coefs.chk to be a vector of length 14 (which seems to line up,
names-wise, with coefs) so 14 seems to be the "right answer".  Where 
does 3 come in?

Previously adding the argument

     "control=glmerControl(check.nobs.vs.nRE="ignore")"

to the call to update() (on Ben's advice) fixed the problems that I was 
having.  No longer is this the case.

What has changed?  Is there any way I can get this to work?

R.

-- 
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276