Package for GLMM with correlation matrix

Hi,

I'm looking for a R package for fitting a generalized linear mixed model
    g(E[y])=X?+Zu
where g is a link function, ? is a p vector of fixed effects, u is a
vector of random effects, X is design matrix, Z is an identity matrix.
In our data, each individual has only one observation for response y.
Because individuals may be correlated in some way, leading to similar
reponses, random effect u is employed to correct for individual
background effect. The variance of u is assumed to be Var[u]=G*?^2, in
which G is a correlation matrix.

glmmPQL in MASS has an argument for correlation. However, I
encountered an error of invalid formula when using a n*n dimension
matrix G (n is the number of individuals, and also the length of y)
for argument 'correlation' in glmmPQL. For example:

y=sample(c(1,0),48,replace=T);x=sample(1:4,48,replace=T);id=1:48;covMat=matrix(rnorm(48*48),nrow=48)
glmmPQL(y~x,random=~1|id,family=binomial,correlation=covMat)

iteration 1
Error in formula.default(object) : invalid formula

The R help documentation for glmmPQL is very compact. No detailed
explanation. Does anybody know how to use correlation matrix in
glmmPQL? Is it good to use glmmPQL for fitting my model?

Thanks.

Alex