Skip to content

extracting p-values from lmer()

3 messages · toka tokas, Renaud Lancelot, Martin Maechler

toka tokas
# 
Dear R users,

I've been struggling with the following problem: I want to extract the Wald p-value
from an lmer() fit, i.e., consider

library(lme4)
n <- 120
x1 <- runif(n, -4, 4)
x2 <- sample(0:1, n, TRUE)
z <- rnorm(n)
id <- 1:n
N <- sample(20:200, n, TRUE)
y <- rbinom(n, N, plogis(0.1 + 0.2 * x1 - 0.5 * x2 + 1.5 * z))

m1 <- lmer(cbind(y, N - y) ~ x1 + x2 + (1 | id), family = binomial, method = "AGQ")
m1


how to extract the p-value for 'x2' from object m1?

Thanks in advance for any hint,
tokas




		
__________________________________________ 

Just $16.99/mo. or less. 
dsl.yahoo.com
Renaud Lancelot
# 
For example:
Generalized linear mixed model fit using AGQ
Formula: cbind(y, N - y) ~ x1 + x2 + (1 | id)
 Family: binomial(logit link)
      AIC      BIC    logLik deviance
 1137.308 1151.246 -563.6541 1127.308
Random effects:
     Groups        Name    Variance    Std.Dev.
         id (Intercept)      3.3363      1.8266
# of obs: 120, groups: id, 120

Estimated scale (compare to 1)  0.8602048

Fixed effects:
              Estimate Std. Error z value  Pr(>|z|)
(Intercept)  0.3596720  0.0070236  51.209 < 2.2e-16 ***
x1           0.2941068  0.0023714 124.025 < 2.2e-16 ***
x2          -0.9272545  0.0100877 -91.919 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
b          se         z P
(Intercept)  0.3596720 0.007023556  51.20939 0
x1           0.2941068 0.002371353 124.02487 0
x2          -0.9272545 0.010087717 -91.91917 0

You might also use the function wald.test in package aod:
Package aod, version 1.1-8
Wald test:
----------

Chi-squared test:
X2 = 15382.2, df = 1, P(> X2) = 0.0

But it is safer to use a likelihood ratio test instead of a Wald test:
Warning message:
IRLS iterations for PQL did not converge
Data:
Models:
m2: cbind(y, N - y) ~ x2 + (1 | id)
m1: cbind(y, N - y) ~ x1 + x2 + (1 | id)
   Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)
m2  4 1149.50 1160.65 -570.75
m1  5 1137.31 1151.25 -563.65 14.192      1  0.0001651 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Best,

Renaud


2005/12/5, toka tokas <tokkass at yahoo.com>:
--
Renaud LANCELOT
D??partement Elevage et M??decine V??t??rinaire (EMVT) du CIRAD
Directeur adjoint charg?? des affaires scientifiques

CIRAD, Animal Production and Veterinary Medicine Department
Deputy director for scientific affairs

Campus international de Baillarguet
TA 30 / B (B??t. B, Bur. 214)
34398 Montpellier Cedex 5 - France
T??l   +33 (0)4 67 59 37 17
Secr. +33 (0)4 67 59 39 04
Fax   +33 (0)4 67 59 37 95
Martin Maechler
# 
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