[R-sig-ME] Another predict/cloglog peculiarity --- lme4 this time.
r@turner @end|ng |rom @uck|@nd@@c@nz
Sun Apr 5 06:07:15 CEST 2020
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
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.
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276
More information about the R-sig-mixed-models