Rejecting Null Too Frequently in Poisson GLMM (lme4, glmmPQL, glmmadmb)

Question

While trying to determine power for a Poisson GLMM, I started by checking the probability of rejecting the null for a given parameter when the null is true (parameter is zero). I kept coming up with a rejection rate of approximately $0.1$ where I expected $\alpha = 0.05$. To check if my programming was faulty, I did the same with a binomial GLMM and a Gaussian LMM. Both of these, however, returned the expected percentage of rejections ($5\%$). Thus, I don't think this is a post for Stack Overflow, but maybe others will disagree.

I figured I'd post the code here and see if anyone can tell me why I'm seeing this unexpected result.

First, a function to simulate data with a single independent variable $X$ and the dependent variable $Y$. The data simulate $100$ individuals with $3$ measurements each. $Y$ is defined as a function of the intercept, $bX$, and the group specific intercept ($g$). Though, $b$ is zero, so $X$ doesn't come into play. The function also fits a GLMM to the data and checks/returns whether $b=0$ should be rejected.

simPow.Pois <- function(j=100, i=3, alpha=.05, b=0, tau=1, refRate=.2) {
    # refRate is referent group rate (intercept)

    # g is group level intercept
    g <- rep(rnorm(j, 0, tau), each=i)

    # group identifies the groupings of the individuals
    group <- rep(1:j, each=3)

    # randomly drawn x
    x <- round(runif(i*j,-.5,5.499)) # approx uniform discrete 0-5

    # DV: a function of intercept, b*x, and group intercept
    y <- exp(refRate + b*x + g)
    y <- rpois(i*j, y)

    # fit the model with one of three options
    #ans <- glmer(y ~ x + (1|group), family=poisson)
    ans <- glmmPQL(y ~ x, random=~1|group, family=poisson)
    #ans <- glmmadmb(y ~ x + (1|group), family = "poisson",link = "log")

    # extract z as [fixed effect] / [SE]
    z <- abs(fixef(ans) / sqrt(diag(vcov(ans))))

    # check if z is too large to believe the null is true
    if(2*(1-pnorm(z['x'])) < alpha) {
        return(1)
    } else { 
        return(0)
    }
}

The function above returns a 1 if the null hypothesis is rejected, and a 0 otherwise.

I then run the following script to do this many times. Unfortunately, this takes a few minutes to run, so I've limited it to 200 iterations. You can expand that if you doubt the result.

require(MASS)
res.pois <- replicate(200,simPow.Pois(b=0))
mean(res.pois)
#[1] 0.105

I appreciate any thoughts on this. Thank you.

Answer 1

I played around with your code and found that it gets closer to 5% when the between individual variance $\tau$ gets smaller. I tried it with $\tau=1, .5, .25, .1$, and got the following results from 1000 replicates:

$\tau$ = 1 --> .103

$\tau$ = .5 --> .07

$\tau$ = .25 --> .051

$\tau$ = .1 --> .051

A possible explanation for this is that Generalized Linear Mixed Models fitted via approximate maximum likelihood based on Penalized Quasi Likelihood (PQL) suffer from biased point estimation when the between individual heterogeneity gets larger because the approximation to the likelihood, based on the Laplace Approximation, is not as accurate in these cases. Have a look at the papers on GLMM by Breslow and Clayton (1993) and by Lin and Breslow (1996) in JASA where this is discussed.

To improve this issue you may consider trying your simulation using an approximation to the likelihood based on adaptive gaussian quadrature (AGQ), as it is known to reduce the bias due to its improved accuracy in approximation to the likelihood. Pinheiro and Bates (1995) in the Journal of Computational and Graphical Statistics discuss comparisons of PQL and AGQ as well as a few other approaches for the nonlinear mixed effects model for a continuous outcome.