How are the p-values in the glmer (lme4) output calculated?

Question

The output of models built using the glmer function in the package lme4 automatically includes p-values. What method does lme4 use to calculate these p-values? I can't seem to find this information anywhere.

Edit: How does lme4 calculate denominator degrees of freedom for GLMMs?

Answer 1

Look at the code of lme4:::summary.merMod. The relevant part is:

p <- length(coefs <- fixef(object))
coefs <- cbind(Estimate = coefs, `Std. Error` = sqrt(diag(vcov(object, 
    use.hessian = use.hessian))))
if (p > 0) {
    coefs <- cbind(coefs, (cf3 <- coefs[, 1]/coefs[, 2]), 
        deparse.level = 0)
    colnames(coefs)[3] <- paste(if (useSc) 
        "t"
    else "z", "value")
    if (isGLMM(object)) 
        coefs <- cbind(coefs, `Pr(>|z|)` = 2 * pnorm(abs(cf3), 
            lower.tail = FALSE))
}

As you see a simple normal approximation is used for GLMMs.

Answer 2

To complete Roland's answer, the pvalues page of lme4 states that:

For glmer models, the summary output provides p-values based on asymptotic Wald tests (P); while this is standard practice for generalized linear models, these tests make assumptions both about the shape of the log-likelihood surface and about the accuracy of a chi-squared approximation to differences in log-likelihoods.

Ben Bolker's excellent GLMM FAQ further states that:

lme4 and similar packages display the Wald Z-statistics for each parameter in the model summary. These have one big advantage: they’re convenient to compute. However, they are asymptotic approximations, assuming both that (1) the sampling distributions of the parameters are multivariate normal (or equivalently that the log-likelihood surface is quadratic) and that (2) the sampling distribution of the log-likelihood is (proportional to) $χ$². The second approximation is discussed further under “Degrees of freedom”. The first assumption usually requires an even greater leap of faith, and is known to cause problems in some contexts (for binomial models failures of this assumption are called the Hauck-Donner effect), especially with extreme-valued parameters.

My understanding is that, since these approximations are "asymptotic", they're particularly inappropriate when the sample size is small (notably when the number of groups in the grouping factor is < 50).