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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
