[R-sig-ME] Another predict/cloglog peculiarity --- lme4 this time.

[Rolf Turner]

On 3/04/20 8:00 pm, Thierry Onkelinx wrote:

> Dear Rolf,
> 
> This is due to the precision of floating point numbers. Let's have a 
> look at what happens when we get the extreme negative difference.
> 
>  > i <- which.min(g(predResp) - predLink)
>  > c(predResp[i], 1 - exp(-exp(predLink[i])))
> 76 76
>   1  1
>  > c(1 - predResp[i], exp(-exp(predLink[i])))
>            76           76
> 2.220446e-16 0.000000e+00
>  > c(-log(1 - predResp[i]), exp(predLink[i]))
>          76         76
>    36.04365 3602.72298
>  > c(log(-log(1 - predResp[i])), predLink[i])
>        76       76
> 3.584731 8.189445

Thanks Thierry.  I think it is at last becoming clear to me.  Sorry for 
taking so long to reply, but I've been thrashing around, trying to 
express things in my own words, in a way that I can understand.

My attempt to do so follows.  Most of you will probably find that what I 
have to say simply obfuscates the issue. You are probably best advised 
to just read what Thierry has written above and ignore my ramblings. 
But just in case anyone is interested, here are my thoughts:

Basically predResp is exactly ginv(predLink); no problems with numerical 
delicacy there.

So when I looked at

     g(predResp) - predLink

I was actually looking at

     g(ginv(predLink)) - predLink

The problem is that g(ginv(x)) is NOT equal to x, when x gets bigger 
than about 3.5, due to numerical delicacy.  (For x < 3.5, g(ginv(x)) is
pretty close to x just as the young and naïve, such as my very good self 
--- :-) --- would expect.)

For x > 3.5, ginv(x) becomes numerically indistinguishable from 1, 
whence g(ginv(x)) is g(1) = Inf.  Or it *would* be except for the fact 
that ginv() (for the cloglog link) is designed so that its values are 
capped at 1 - .Machine$double.eps = 2.220446e-16 (on my machine at least).

Consequently g(ginv(x)) is capped at g(1 - 2.220446e-16) = 3.584731.
So as x grows it quickly outstrips g(ginv(x)).

When x is close to 3.5 g(ginv(x)) is simply a bit numerically wobbly, 
and not necessarily smaller than x. (Which is why I got a negative 
value, explicitly -0.00281936, for the lower bound of the range of
g(predResp) - predLink.

Another person, who replied to me off-list, suggested plotting
g(predResp) ~ predLink .  This is indeed illuminating.

I think that the main thing to take away from all this is that numerical 
delicacy is always lurking in the bushes, and care should always be 
taken, especially when logs and exponentials are involved.  It is 
perilous to assume that *algebraic* identities, like g(ginv(x)) = x,
will turn out to be *numerically* true.

Thanks again to Thierry.

cheers,

Rolf

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