Skip to content

Cross validated log likelihood, redux.

3 messages · Rolf Turner, Ben Bolker

Original
Rolf Turner
# 
My apologies for continuing to pester the list with questions about this 
issue, but I urgently need an answer to my most recent question.

That question was contained in a postscript to a reply to D. Rizopoulos, 
and thereby may have been overlooked.  Consequently I am re-posting this 
question.  (Again I apologise for taking up bandwidth.)

In summary, the situation is as follows:

* I am trying to calculate the log likelihood of a "new" data set on the 
basis of a model fitted to a different data set.

* Ben Bolker showed me how to do this using glmer() (from the lme4 
package) and after some to-ing and fro-ing I got his recipe to work.

* Dimitris Rizopoulis also showed me how to do this using the 
mixed_model() function from the GLMMadaptive package.  Again, after much 
delay and after more to-ing and fro-ing, I got Prof. Rizopoulis's advice 
to work.

* I then wanted to cross-check the value obtained from mixed_model() 
with that obtained from glmer(), so I re-ran the glmer() based code.
Lo and behold, that code threw an error (where it had not before done so).

I would really like to be able to use both methods (i.e. that based on 
mixed_model() and that based on glmer()).  So I would like to figure out 
what is going wrong --- or what I am doing wrong --- in the case of the 
glmer() approach.

I have attached a source-able script demo.glmer.txt to demonstrate what 
happens, and have also attached the (simulated) data set X.txt which the 
script uses.

If you place demo.glmer.txt and X.txt in your working directory and
source("demo.glmer.txt") you will get the following message and error 
message:
As I said in my (possibly overlooked) postscript in a previous posting 
on this issue:

This code *worked* previously!!!

I thought that the dropped coefficient might be the problem, so in the 
demo script 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 as produced by the script) so 14 seems to be the 
"right answer".  Where does 3 come in?

Initially (way back when) I got an error *something* like this, but then
Ben Bolker advised me to add the argument

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

to the call to update() and that fixed the problem that I was having. 
But the fix no longer seems to be "operative". :-)

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

Thanks for any pearls of wisdom.

cheers,

Rolf Turner
Ben Bolker
# 
Hmmm. The proximal problem is that the updated deviance function
thinks it needs coefficients for all of the parameters (theta = 3 +
beta=12), not just the theta parameters.

  It's not obvious to me (I know this might be code I wrote in the first
place!) why you're trying to pass the full theta+beta coefficient vector
to a model defined with nAGQ=0 (which estimates the beta coefficients as
part of the penalized weighted resid sum of squares (pwrss)
computation).  Can you explain/remind us of the reason for the mismatch
here?

  cheers
    Ben Bolker
On 2019-08-16 7:55 p.m., Rolf Turner wrote:
Rolf Turner
# 
Hi Ben.  Thanks for getting back to me on this.
On 17/08/19 1:43 PM, Ben Bolker wrote:
The short answer is, I haven't got a clue.  (Or perhaps, to put it 
differently, that I'm clueless. :-) )

Basically I have just been (mindlessly) following instructions from your 
very good self.  Initially (20 April 2019) you said:
I tried that, and it didn't work.  (Errors were thrown.) Then you 
updated the recipe (on 26 April 2019):
having remarked that there were *two* issues.  (One issue involved the 
specification of "coefs", the other involved having to "override some 
checking that glmer does".)

In my implementation of the code to do a cross-check with the results 
from mixed_model() I neglected to take cognizance of the first issue. 
I.e. I left the assignment of "coefs" as:
whereas I should have changed it to
i.e. leaving out the "beta" term.  Duhhhh!!!

Having corrected my error I find that the code now works as before.  Not 
surprisingly.

Thank you for asking for an explanation of what I was doing, which 
prompted me to go back over the story with the result that I spotted the 
botch-up that I had made.

(It's very difficult to spot errors when you are flying intellectually 
blind.)

Bottom line:  the harmony of the universe now seems to have been 
essentially restored. :-)

Thanks again.

cheers,

Rolf

P. S.  In case anyone is interested, the log likelihood of the 
"validation" data set that I get from mixed_model() in the example that 
I provided is -35.26048.

That from glmer() is -21.37535.

I guess that one can say that they are "in the same ball park".

Since I have no real understanding of the underlying intricacies, I have 
no idea what the implications of the difference between the two answers 
might be.

R. T.