Is it possible to fit a crossed random effects model using glmmPQL and, e.g. random=pdBlocked(list(pdIdent(~a-1),pdIdent(~b-1))) as the random part of the model? I seem to be getting errors on R using glmmPQL in this way, although if I try something very similar using 'lme' instead of 'glmmPQL', it works fine. Does the data have to be of a groupedData form for it to work (as it has to be for lme) because glmmPQL doesn't seem to work with data of the form groupedData?