# Re: [R] extracting p-values from lmer()

From: Renaud Lancelot <renaud.lancelot_at_gmail.com>
Date: Tue 06 Dec 2005 - 18:09:35 EST

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@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@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
>

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

______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help