[R] extracting p-values from lmer()

[Renaud Lancelot]

For example:

> m1
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
>
> 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)
                     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:

> library(aod)
Package aod, version 1.1-8
> wald.test(Sigma = vc, b = b, Terms = 2)
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:

> # LRT to test the coef associated with x1
> m2 <- lmer(cbind(y, N - y) ~ x2 + (1 | id), family = binomial, method = "AGQ")
Warning message:
IRLS iterations for PQL did not converge
> anova(m1, m2)
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>:
> 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
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>


--
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