[R] pdIdnot / logLik in glmmPQL

[ptr2003@columbia.edu]

Dear R users,

I have been using the pdMat class "pdIdnot" (from the mgcv
package)instead of "pdIdent" to avoid overflow in GLMM fits with
the MASS package function glmmPQL, of the following form:

fit1 <- glmmPQL(fixed=y0~-1+xx0, random=list(gp=pdIdent(~-1+zz0)),
                      family=binomial) # vulnerable to overflow
fit2 <- glmmPQL(fixed=y0~-1+xx0, random=list(gp=pdIdnot(~-1+zz0)),
                      family=binomial) # overflow-proof

In instances in which fit1 does *not* lead to overflow, the result
sometimes differs from fit2.  This leads me to two questions.

1. Does anyone have any thoughts on what might cause such a
discrepancy?

2. Given two discrepant fits, I would like a way to choose the
better one.  If my reading of Breslow and Clayton's 1993 paper
(specifically, their equation 12) is correct, at convergence, the
profile quasilikelihood should be approximately equal to the log
likelihood from the last linear mixed model fit by the algorithm. 
If so, the $logLik component of the lme object produced by glmmPQL
should approximate the quasilikelihood I am trying to maximize.  In
short: according to this argument, the glmmPQL fit with higher
"logLik" should be the better one.

And yet, some previous postings seem to indicate that the "logLik"
of an object produced by glmmPQL cannot be interpreted in terms of
likelihood, quasi or otherwise.  If the above made sense (or even
if not), is there anyone who could kindly speak to this point?

Any help with either of the above questions would be greatly
appreciated.

Phil Reiss
PhD candidate, Dept. of Biostatistics
Columbia University, New York
ptr2003 at columbia.edu