[R] extracting p-values from lmer()
Martin Maechler
maechler at stat.math.ethz.ch
Tue Dec 6 10:03:31 CET 2005
>>>>> "Renaud" == Renaud Lancelot <renaud.lancelot at gmail.com>
>>>>> on Tue, 6 Dec 2005 08:09:35 +0100 writes:
Renaud> For example:
....
>> vc <- vcov(m1, useScale = FALSE)
>> b <- fixef(m1)
>> se <- sqrt(diag(vc))
>> z <- b / sqrt(diag(vc))
>> P <- 2 * (1 - pnorm(abs(z)))
>>
>> cbind(b, se, z, P)
Renaud> b se z P
Renaud> (Intercept) 0.3596720 0.007023556 51.20939 0
Renaud> x1 0.2941068 0.002371353 124.02487 0
Renaud> x2 -0.9272545 0.010087717 -91.91917 0
I still see much too many uses of "1 - p<dist>(...)"
which in cases as the above case leads to complete loss of
accuracy (1 - 1 = 0) -- well actually the above case is too
extreme to make any difference; but let me explain the general principle:
Though the loss is usually no problem for decision making based on
P-values, it is unnecessary:
One of the (extra) features of R are the arguments 'lower.tail'
and 'log.p' of all the p<dist>() functions -- which (in not yet
quite all cases) allow avoid precision loss.
E.g.,
> 1 - pnorm(c( 6,8,10,20))
[1] 9.865877e-10 6.661338e-16 0.000000e+00 0.000000e+00
> pnorm(c(6,8, 10,20), lower.tail=FALSE)
[1] 9.865876e-10 6.220961e-16 7.619853e-24 2.753624e-89
BTW, example(pnorm) ends in two plots which show the
advantage of using 'log.p' for additional precision gain
e.g. for log-likelihood computation.
Martin Maechler, ETH Zurich
More information about the R-help
mailing list