Hello R-Users, I am using glmmPQL (library MASS) on time-series count data that are aggregated by zipcode. My model includes natural cubic splines for season and day-of-week as fixed effects and random intercept term for zipcode. I need to extract the fitted values AND the standard error of the fit to compute confidence intervals for the fitted value. I am able to extract the fitted, but the predict function doesn't work like it does for GLM [predict(glm_obj, se.fit=T)]. Is there an equivalent for glmmPQL? Or is there some other way I can extract se? Would also be good to know if this is not possible and there is an alternate glmm model I can consider. Thank you! Ramona Ramona Lall, PhD City Research Scientist Bureau of Communicable Diseases New York City Department of Health and Mental Hygiene [[alternative HTML version deleted]]