[R] obtaining means/SD after fitting a mixed model

[Pablo Inchausti]

Hello,
I am running (non-balanced) mixed models (using library lme4) such as :
model1<-glmmPQL(Y~grouping variable+ covariate, random=~1|yr/month, 
data=dat, family= gaussian or poisson)
and besides the usual output in terms of the statistical significance of 
the fitted coefficients, I'd like to obtain the adjusted means and standard 
deviations (as one obtains after fitting a lm model). My problem is that I 
cannot find a command equivalent to "model.table (model, grouping variable, 
mean/sd)" that would apply to mixed models. Any suggestion would be 
appreciated.
Thanks in advance

Pablo