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

[Rolf Turner]

Dear Prof. Rizopoulos,

I hope that this gets through to you on your travels and that I am not 
disrupting the even tenor of your ways *too* badly. :-)

I am *finally* getting back to my efforts to calculate the log 
likelihood of a new ("validation") set of data on the basis of a model 
fitted to a different ("training") set of data.

You suggested a procedure for doing this using the mixed_model() 
function from the GLMMadaptive package.  After a bit of running-in time, 
getting to know the mixed_model() function, I (finally!) set about 
trying to follow your advice.

On 29/04/19 4:25 PM, D. Rizopoulos wrote:

> 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)

I am a bit mystified by the "iter_EM = 0" and "iter_qN_outer = 0" in the 
foregoing.  I presume that they really shouldn't be there, but belong 
only in the next step.  But I could be completely out to lunch here, 
since I don't know what I'm doing at all!

I presume that "Ts" should be the "training set" of data.

> # 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)

I presume (yet again!) the "Ts" in the foregoing should be "VS", the new
"validation" set of data.  If not, then I'm *really* not understanding 
anything that is going on here.

What I am conjecturing is that the foregoing code "fits" a model to the
"validation" data without actually doing any fitting.  I.e. it just 
"stays at the starting values" (since it does no iterations) thereby 
getting the same model as that encompassed by "fm".  If this is not what 
the code is doing, can you please explain what it *is* doing?

For "beta" in the initial_values list I used fixef(fm) and for "D" I 
used fm$D.  If this isn't right, then I need more instruction.

Anyhow:  I put all this together in a script ("demo.txt") which is 
attached, along with a demonstration data set X.txt.

If these file are save into your working directory you should be able to
source("demo.txt") and see what I saw:

> Error in if (use.names && nt[i] == nc[i]) dQuote(nt[i]) else i : 
>   missing value where TRUE/FALSE needed
> In addition: Warning messages:
> 1: glm.fit: fitted probabilities numerically 0 or 1 occurred 
> 2: glm.fit: fitted probabilities numerically 0 or 1 occurred 

I presume (am I presumptuous or what?) that the warnings are not too 
much to worry about, but the error flummoxes me.  Clearly I've got 
something wrong, but I have insufficient insight to see what it is.

Can you (or possibly someone else) enlighten me?

Thanks.

cheers,

Rolf

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

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: demo.txt
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20190814/ad202abb/attachment.txt>

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: X.txt
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20190814/ad202abb/attachment-0001.txt>