[R-sig-ME] lmer vs glmmPQL

[Roger Levy]

Ben Bolker wrote:
> Daniel Ezra Johnson wrote:
>> Regarding this thread, what about the method of fitting nested models
>> and using anova() to estimate a p-value.
>>
>> For example:
>>
>> mod2 = lmer(y ~ genotype + (1|block), family = binomial, data)
>> mod0 = lmer(y ~ (1|block), family = binomial, data)
>>
>> anova(mod0,mod2)
>>
>> How does that p-value (which is one value for the whole term
>> "genotype") relate to the individual coefficient p-values derived from
>> the Z-scores inside summary(mod2)?
>>
>> Thanks,
>> Dan
> 
>    The Z-scores are Wald tests.  The anova results are likelihood
> ratio tests.  The latter are in general more reliable but are known
> to be *un*reliable for small-sample (i.e.
> small-number-of-random-effects-levels) LMMs (Pinheiro and Bates).
> Wald tests are not *known* to be unreliable in the small-sample
> case, but I believe that is a statement of ignorance rather than
> a statement that they're OK ...

This is an interesting issue -- as Ben points out, likelihood-ratio 
tests are known to be anti-conservative for tests between linear mixed 
models differing only in fixed-effects structure.  I think the 
best-known demonstration of this can be found in Pinheiro & Bates, 2000, 
pages 87-92.  F-tests are less problematic.

For logit mixed models, presumably the results will depend a bit on how 
fitting is done because the likelihood has to be approximated.  For the 
Laplace approximation currently in use in lme4, I've done some 
simulations and generally found that likelihood-ratio tests are still 
anti-conservative; tests based on Wald Z are not as anti-conservative, 
though they don't seem quite nominal.  You can do simulations along 
these lines:

set.seed(2)
invlogit <- function(x) {
   temp <- exp(x)
   return(temp/(1+temp))
}
L <- 10
REP <- 4
N <- 4*10
fac <- gl(L,REP)
f <- function(ignore) {
   x <- runif(N)
   r <- rnorm(L)
   y <- rbinom(N,1,invlogit(r[fac]))
   m0 <- lmer(y ~ (1 | fac), family="binomial")
   m1 <- lmer(y ~ x + (1 | fac), family="binomial")
   temp <- pnorm(fixef(m1)[2]/sqrt(vcov(m1)[2,2]))
   p.z <- 2*min(temp, 1-temp)
   p.lr <- 1-pchisq(as.numeric(2*(logLik(m1)-logLik(m0))),1)
   return(c(p.z,p.lr))
}
K <- 5000
p.vals <- sapply(1:K,f)
apply(p.vals, 1, function(x) sum(x<0.05))/K
# this gave me 0.0562 for Wald Z, 0.0608 for LRT

# Is Wald Z anti-conservative?  The below says "yes" with p=0.044
2*(1-pbinom(sum(p.vals[1,]<0.05), K, 0.05))


Unlike what Ben states above, however, it's my intuition (and informal 
observation on the basis of simulations) that what determines how 
anticonservative the LRT is is not the number of random effect levels 
but rather the ratio between the degrees of freedom involved in the LRT 
to the residual degrees of freedom in the more general model.  The more 
residual degrees of freedom (basically, the more observations), the less 
anti-conservative the LRT is.  This can be assessed for a specific 
dataset by fitting a null-hypothesis model, repeatedly generating new, 
artificial responses from it, and looking at the distribution of LRTs on 
the artificial responses (along the lines of P&B 2000).

Roger