How do I plot predictions from new data fit with g

I have fit my discrete count data using a variety of functions for comparison. I fit a GEE model using geepack, a linear mixed effect model on the log(count) using lme (nlme), a GLMM using glmer (lme4), and a GAMM using gamm4 (gamm4) in R.

I am interested in comparing these models and would like to plot the expected (predicted) values for a new set of data (predictor variables). My goal is to compare the predicted effects for each model under particular conditions (x variables). Of particular interest is the comparison between marginal (GEE) and conditional estimates.

I think my main problem might be getting the new data in the correct form with the correct labels and attributes and such. I am still very much an R novice and struggle with this stuff (no course on this at my university unfortunately).

I currently have fitted models

gee1 lme1 lmer1 gamm1

and can extract their fixed effect coefficients and standard errors without a problem. I also don't have a problem converting them from the log scale or estimating confidence intervals accounting for the random effects.

I also have my new dataframe newdat which has 365 observations of 23 variables (average environmental data for each day of the year).

I am stuck on how to predict new count estimates from this. I played around with the model.matrix function but couldn't get it to work. For example, I tried:

mm = model.matrix(terms(glmm1), newdat) # Error in model.frame.default(object, 
                                        # data, xlev = xlev) : object is not a matrix
newdat$pcount = mm %*% fixef(glmm1)

Any suggestions or good references would be greatly appreciated. Can anyone help with the error above?

Getting predictions for lme() and lmer() is documented on http://glmm.wikidot.com/faq