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
> 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 ______________________________________________ R-help@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
Received on Tue Dec 06 18:14:53 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:41:29 EST