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 <[EMAIL PROTECTED]>:
> 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
>
> ______________________________________________
> [email protected] 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
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
