[R-sig-ME] extracting fitted means from nlme and glmer

[Jenny Hutchinson]

Hello R mixed-modelers,

I am a mixed models and R novice and probably know just enough to be
dangerous, if not entirely confusing, so please forgive any misuse of
terminology or incorrect posting etiquette.  I've searched the
archives and looked elsewhere online and cannot find a solution to my
problem.  Here's what I'm trying to do.

I am interested in extracting fitted means of my response variables
'pollen grain count' and fruit weight' that take into account
variability in covariates from an lme model (nlme package).  I will
also be trying to do the same for my binary response variable
'proportion fruit set' from a glmer model (lme4 package).  I have five
covariates (height, width, density, distance to tree, and distance to
non-coffee) and three factors with three levels each (canopy cover,
farm, and treatment).  My model for grain count is as follows:

model1=lme(sqrt(GC)~Farm*TRT+Height+Width+Density+DT+DNC+CC,data=gc,random=~1|Bush,na.action=na.omit)

Treatments were nested within bushes so 'Bush' is my random effect and
the response was square root transformed.  I want to extract and graph
(using the plotrix package) the fitted means of 'pollen grain count'
with their appropriate SE's, realizing that after back-transformation
the intervals will be asymmetric.  My statistics advisor wrote a
little code that works for a simple data set he has with one factor
(species) and one covariate (length) that predicts 'eggs'.  This is
not working when I try to conform it to my data and I suspect that it
is due to the inclusion of the random effect 'Bush'. Here is my code
so far (again apologies for inelegance):

gc=read.csv("graincount.csv")
head(gc) #Grain count is response; also contains 5 covariates and 3 categorical
attach(gc)
model1=lme(sqrt(GC)~Farm*TRT+Height+Width+Density+DT+DNC+CC,data=gc,random=~1|Bush,na.action=na.omit)
#Fit model
summary(model1)
Covmeanheight=mean(Height, na.rm=TRUE) #mean of covariates
Covmeanwidth=mean(Width, na.rm=TRUE)
Covmeandensity=mean(Density)
CovmeanDT=mean(DT, na.rm=TRUE)
CovmeanDNC=mean(DNC, na.rm=TRUE)
covcolumnheight=rep(Covmeanheight,27) #column vector containing
covariate mean (27 times because there are 27 combinations of groups)
covcolumnwidth=rep(Covmeanwidth,27)
covcolumndensity=rep(Covmeandensity,27)
covcolumnDT=rep(CovmeanDT,27)
covcolumnDNC=rep(CovmeanDNC,27)
TRTcolumn=c("Open","Bagged","Crossed","Open","Bagged","Crossed","Open","Bagged","Crossed","Open","Bagged","Crossed","Open","Bagged","Crossed","Open","Bagged","Crossed","Open","Bagged","Crossed","Open","Bagged","Crossed","Open","Bagged","Crossed")
#Column vector containing 3*9 group destinations
CCcolumn=c("Sun","Sun","Sun","Shade","Shade","Shade","Intermediate","Intermediate","Intermediate","Sun","Sun","Sun","Shade","Shade","Shade","Intermediate","Intermediate","Intermediate","Sun","Sun","Sun","Shade","Shade","Shade","Intermediate","Intermediate","Intermediate")
Farmcolumn=c("KP09","KP09","KP09","KP09","KP09","KP09","KP09","KP09","KP09","KP10","KP10","KP10","KP10","KP10","KP10","KP10","KP10","KP10","BH10","BH10","BH10","BH10","BH10","BH10","BH10","BH10","BH10")
newdata=data.frame(covcolumnheight,covcolumnwidth,covcolumndensity,covcolumnDT,covcolumnDNC,TRTcolumn,CCcolumn,Farmcolumn)
names(newdata)= names(gc)[2:9] #column one contains 'Bush' ID
newdata
fittedmeans=predict(model1,newdata,se.fit=TRUE,type="response")

When I get to this point I get the error message:

Error in predict.lme(model1, newdata, se.fit = TRUE, type =
"response") : Cannot evaluate groups for desired levels on "newdata"

His code is as follows so you can see how I've converted it to fit my data set:

f=read.csv("fecundity.csv")
head(f) #Eggs is response; length (covariate) and species are predictors
attach(f)
Length=log(Length) #I had to log-transform my x variable
model=lm(log(Eggs)~Length*Species) #Fit model
summary(model)
Covmean=mean(Length) #Mean of (transformed) covariate
covcolumn=rep(Covmean,2) #Column vector containing the covariate mean
(twice because there are two groups))
Spcolumn=c("YCT","RHT") #Column vector containing the two group designations
newdata=data.frame(covcolumn,Spcolumn) #Combine the two columns into a dataframe
names(newdata)= names(f)[1:2]  #Put correct predictor names
newdata  #This is now a data frame with the covariate mean as "Length"
and the two species in "Species"
fittedmeans=predict(model,newdata,se.fit=TRUE,type="response")
#Predict mean response of the groups at this mean covariate value
fittedmeans$fit  #Fitted means
fittedmeans$se.fit  #SE of fitted means
intervals=data.frame(fittedmeans$fit-2*fittedmeans$se.fit,fittedmeans$fit,fittedmeans$fit+2*fittedmeans$se.fit)
names(intervals)=c("mean-2SE","mean","mean+2SE")
intervals  #Mean +/- 2 SE in transformed units
backints=exp(intervals) #Backtransform, given the original log
transform in this problem
plotCI(c(0,1),backints[,2],li=backints[,1],ui=backints[,3],ylab="Eggs")
#plots back-trasnformed fitted means +/- 2SE
#But I don't know how to get the group labels on the x-axis

What am I doing wrong?  Is there an easier method than this?  I would
like to post my data so the error can be replicated but am unsure how
to do so to the list. Sorry for the lengthy posting!

Any help or advice would be greatly appreciated!  I am already greatly
indebted to this list as it has already helped solve many problems up
until now!


Jenny Lynn

--
Jenny Lynn Hutchinson
M.A. Candidate
Biological Sciences
Humboldt State University
707-407-6244
jennylynnhutch at gmail.com