Another predict/cloglog peculiarity --- lme4 this time (PS)

cllinv <- function(eta) <- 1-exp(-exp(eta))

options(digits=17)

rlinkinv <- make.link('cloglog')$linkinv
cllinv <- function(eta)1-exp(-exp(eta))

cll <- function(mu)log(1-log(mu))

mat <- matrix(nrow=6, ncol=3)
colnames(mat)<- c('1-exp(1-exp(eta))','R linkinv()', 'Calc as 1 -')
rownames(mat) <- paste(seq(from=3.5, to=4, by=0.1))
i <- 0
for(eta in seq(from=3.5, to=4, by=0.1)){
  i <- i+1
  mat[i,] <- c(cllinv(eta), rlinkinv(eta), exp(1-exp(eta)))
}
print(mat)

cll(rlinkinv(3.7))

On 5/04/2020, at 16:07, Rolf Turner <r.turner at auckland.ac.nz> wrote:


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