0

Hi I recently run into the problem of GLMM estimation, specifically, the bias of the variance components (and the fixed-effects) estimation using the PQL method. Two key references are Brewlow&Clayton 1993 and Lin&Breslow 1995, which both claims that PQL under-estimates the variance of the random-effect, see Section 4.3 and Table 1 of the first paper, and Eq (3) and Fig 4 of the second paper. Both papers are highly cited so this should be well accepted.

However when I tried to replicate the simulation, the results are totally different. I coded in R and used MASS::glmmPQL(). Below are some quick code.

## Brewlow&Clayton 1993, Table 1, D=D1, m=1
require(MASS)
require(nlme)
require(lme4)

set.seed(123)
sim.breslow <- matrix(NA, 200, 4+1)
expit <- function(a) 1/(1+exp(-a))    
K = 100         # 100 clusters
nn = rep(7, K)  # size of each cluster
N = sum(nn)
pid <- rep(1:K, nn)  # cluster id
xx = rep(c(1,0), each=K*7/2)
tt = rep(seq(-3,3), K)
X = cbind(tt, xx, tt*xx)
for(ii in 1:200){     # 200 reps
  g.t <- rnorm(K, 0, 1)   # random intercept
  lp = -2.5 + X %*% c(1,-1, -0.5) + rep(g.t, each=7)  # linear predictor
  y = rbinom(N, 1, expit(lp))
  dd = data.frame(y,X,pid=pid)
  names(dd)[2:4] <- c('tt', 'xx', 'xt')
  fit.pql = tryCatch(glmmPQL(fixed = y ~ tt+xx+xt, random = ~1|pid, family = 'binomial', data = dd, verbose=F),
                     error = function(e) NA)
  if(length(fit.pql)!=1){
    b.pql <- fit.pql$coef$fixed
    D.pql <- as.numeric(VarCorr(fit.pql)[1,1])
  }else{
    b.pql <- rep(NA, 4)
    D.pql <- rep(NA, 1)
  }
  sim.breslow[ii,] <- c(b.pql, D.pql)
cat('.')
}

res = apply(sim.breslow, 2, mean) 
res = cbind(c(-2.5, 1,-1, -0.5, 1), res)
row.names(res) <- c('a0', 'a1', 'a2', 'a3', 's00')
colnames(res) <- c('true', 'PQL')
round(res, 2)
#     true   PQL
# a0  -2.5 -2.61
# a1   1.0  1.04
# a2  -1.0 -0.93
# a3  -0.5 -0.52
# s00  1.0  1.72

Notice that the var component s00 is over-estimated (1.72 compared to the truth 1.0), this is different from the paper claims (under-estimated as 0.68 while the truth is 1.0). I also tested other settings in the two papers and the same conclusion, I can post more simulation code if necessary. I checked the source code of glmmPQL and seems no problem.

So my questions are

  1. Why glmmPQL result differs from the conclusion the references claimed? Did I miss anything?
  2. Are there any other functions I can use to cross-validate? MASS::glmmPQL is the only one I can find to do PQL.
  3. BTW, anyone aware of any packages that do the bias correction of Lin&Breslow 1995?

Thanks!

Jason Luo
  • 3
  • 1
  • 5
  • The method by `MASS` calls the implementation by `lme`, so potentially some of the deviation is from small differences in algorithm (although, `PQL` is a very well defined method). If you want to be certain implementing the different parts of the algorithm is not a far stretch to do. One thing talking for the slight difference in algorithm, is that if you repeat your test for `m = 1, 2, 4, 8` we see that the variance estimate converges towards the true estimate, however from the opposite direction to B&C's observation. It might even come down to initial choice of `sigma_00` – Oliver Dec 06 '20 at 12:28
  • In their article they choose very specific starting points while I'm somewhat certain that `lme` initializes these uniformly. A thing talking against this thought, is of course the function should converge to a single optimum (as it is convex). But I am just spitballing a few suggestions. – Oliver Dec 06 '20 at 12:29
  • @Oliver thanks very much for your comments (I realized this probably is not a good question on SO). I repeated my test for m=2,4,8 and the estimated `sigma_00` is 1.12, 0.93, 0.93, so basically as you envisioned. However I'm still confused about the `lme` part. As I understand, PQL converts the problem to LMM estimation, which is quite easy to solve via profile likelihood, and not depend on actual implementation (I modified the `MASS` code to use `lme4` instead of `lme` and got almost identical results). – Jason Luo Dec 07 '20 at 23:31
  • It should be close to identical results if you convert the code base to use `lme4`. The two are optimizing the same likelihood function. As for *PQL* it doesn't "convert the problem to LMM estimation" but something quite similar. Using some simplifications we can obtain a marginal likelihood that is (almost) "identical" to the conditional likelihood. Basically it means we're optimizing the fixed effect **and** the random effects using the conditional likelihood. Just like optimizing GLM models, we can use Iterative least squares (eg. LM optim) to find our parameters. – Oliver Dec 08 '20 at 13:11
  • Breslow and Clayton take a similar approach, but maybe there is a difference that is not "simple to spot". You are indeed right it might be better suited for crossvalidation. Honestly I am at a loss why there would be a difference. – Oliver Dec 08 '20 at 13:12
  • @Oliver I checked the SAS version of PQL, the results are different, so probably MASS has a different implementation. Will explore it later. Feel free to ANSWER my question and I will take it. Thanks! – Jason Luo Dec 13 '20 at 04:35

0 Answers0