[R] how to get means and confidence limits after glmmPQL or lmer

[Don Driscoll]

R,

I want to get means and confidence limits on the original scale for 
the treatment effect after running a mixed model.

The data are:

response<-c(16,4,5,8,41,45,10,15,11,3,1,64,41,23,18,16,10,22,2,3)

block<-factor(c("A","A","A","A","I","I","I","I","N","N","N","P","P","P","P","P","P","S","S","S"))

treatment<-factor(c("b10","b10","b20","b20","b10","b1","b1","b20","b10","b20","b1","b10","b10","b1","b1","b20","b20","b10","b20","b1"))



library(MASS)
mm<-glmmPQL(response~ treatment, random = ~ 1 | block,family = poisson)


antab<-anova.lme(mm)
antab

This tells me that response is significant.  All I want to do now is 
get the means and CIs for the levels of treatment on the original 
scale.  Below I describe some of the things I have tried.



nn<-intervals(mm)#to get confidence intervals on estimates
oo<-nn$fixed
oo[2,]<-oo[1,]+oo[2,]
oo[3,]<- oo[1,]+oo[3,]
pp<-exp(oo)

But the means are not the right means (see speciesmeans below), the 
intercept does not correspond with level 1 of treatment (b1), and the 
confidence limits seem too big.  Also, I'm not sure if I've 
backtransformed correctly.

I've also tried directly calculating the means and CIs, but I'm not 
sure if the first line below is actually the right variance to use in 
estimating the CIs.


variance<-mm$sigma  #IS THIS THE RIGHT VARIANCE?
  speciessize<- tapply(species, burnagecat1 ,length)
  speciesmeans<- (tapply(predict(mm, type="response"), burnagecat1 ,mean))
  upperci<- exp(log(speciesmeans) + 
qt(0.95,residdf)*sqrt(variance/speciessize))
lowerci<- exp(log(speciesmeans)- qt(0.95,residdf)*sqrt(variance/speciessize))

result<-cbind(lowerci, speciesmeans, upperci)
result


I also explored lmer, but haven't found a way to extract means and 
CIs on the original scale for the levels of treatment.

library(lme4)
dlm<-lmer(response~ treatment+ (1|block),family=poisson)
summary(dlm)

A suggestion from Spencer Graves in 2006 to use mcmcsamp and 
subsequent manipulations to get a distribution for CI estimation fails:

 > mcmcsamp(dlm, n=1000, verbose)
Error in .local(object, n, verbose, ...) : Update not yet written

I've also tried predict for the glmmPQL model (mm) and for the lmer 
model (dlm):


ses<-predict(mm, 
data.frame(treatment=c('b1','b10','b20')),type="response", se.fit=TRUE)

Error in predict.lme(object, newdata, level = level, na.action = na.action) :
   Cannot evaluate groups for desired levels on "newdata"


ses2<-predict(dlm, 
data.frame(treatment=c('b1','b10','b20')),type="response", se.fit=TRUE)
Error in UseMethod("predict") : no applicable method for "predict"



Any advice on how to get the means and CIs on the original scale 
after fitting my simple glmm using either glmmPQL or lmer will be 
much appreciated.




Don
(using R version 2.9.0)

Dr Don Driscoll
Fenner School of Environment and Society
WK Hancock Building 43
Australian National University
CANBERRA ACT 0200
02 6125 8130